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QUARTERLY OF APPLIED MATHEMATICS 


Vol. IV JULY, 1946 No. 2 


Harry Bateman* 
29 May, 1882—21 January, 1946 


In the sudden death (from coronary thrombosis) of Harry Bateman while en route 
to New York, near Milford, Utah, mathematics in the United States lost its outstand- 
ing representative of the British School of the generation just closing. Like his con- 
temporaries and immediate predecessors among Cambridge mathematicians of the 
first decade of this century, before the new regulations for the Mathematical Tripos 
took effect, Bateman was thoroughly trained in both pure analysis and mathematical 
physics, and retained an equal interest in both throughout his scientific career. In 
bare outline the relevant details of his life are as follows: 

Harry Bateman was born at Manchester, England, 29 May, 1882, a son of Samuel 
and Marnie Elizabeth (Bond) Bateman, and received his secondary education at the 
Manchester Grammar School. Bateman ascribed much of his subsequent success at 
Trinity College, Cambridge, to the excellent instruction he received at the school. 
In 1903 he was (bracketed) Senior Wrangler in the Tripos, and took his B.A. degree, 
proceeding to the M.A. in 1906, having been a Smith’s Prizeman in 1905. From 1905 
to 1911 he was a Fellow of Trinity College: the year 1905-06 was spent in study at 
Géttingen and Paris. From 1906 to 1907 he was a lecturer at the University of Liver- 
pool, and from 1907 to 1910 a reader at the University of Manchester. He came to 
the United States in 1910 (he later became a naturalized U. S. citizen), as a lecturer 
at Bryn Mawr College, where another English mathematician, the late Charlotte 
Angas Scott, was the efficient and scholarly head of the mathematics department. 
In 1912, he went to the Johns Hopkins University as a Johnston scholar for three 
years, incidentally taking his Ph.D. (a curious proceeding for a mathematician of his 
eminence) in 1913. From 1915 to 1917 he was a lecturer at Johns Hopkins, and in 
1917 he accepted the position which he held till his death, a professorship of mathe- 
matics, physics, and aeronautics at the then recently organized California Institute 
of Technology. He was a member of the American Mathematical Society (vice- 
President, 1935, Gibbs lecturer, 1943), the American Physical Society, the American 
Acoustical Society, the American Philosophical Society, the British Association for 
the Advancement of Science, the London Mathematical Society, the National (U. S.) 
Academy of Sciences, and a Fellow of the Royal Society (London). He is survived 
by his wife, Ethel (Horner Dodd) Bateman, and his daughter, Joan; a son died in 
childhood. 


* Professor Bateman was a member of the Board of Collaborators of the Quarterly of Applied Mathe- 
matics from its foundation to his lamented decease. 
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Bateman was an almost unique combination of erudition and creativeness. It is 
most unusual for a mathematician to have the extraordinary range of exact knowledge 
that Bateman had, and not be crushed into sterility by the mere burden of an oppres- 
sive scholarship. But, as his numerous publications testify, Bateman retained his 
creative originality till his death. In pure mathematics, his dominating interest was 
in the analysis that has developed from classical mathematical physics. His technical 
sixill in this broad field was unrivalled. His numerous contributions to mathematical 
physics are marked by a vivid, at times almost romantic, imagination. Students of 
the history of general relativity will find much of interest in some of his papers on 
electromagnetism. 

A singularly modest and gentle man, Bateman was always ready to place his skill 
and his knowledge at the disposal of others, with no thought of personal credit. War 
work absorbed most of his time during the last four years of his life; and it is to be 
regretted that the incessant correspondence in connection with such work prevented 
him from putting the finishing touches to what he regarded as his most useful con- 
tributions to mathematical scholarship: an exhaustive work on definite integrals, and 
a critical census of all the special functions that have been considered in mathematics. 
If these works can be put into shape for publication, they will form a lasting memorial 


to Harry Bateman. 
E. T. BELL 


April, 1946. 


LIST OF PUBLICATIONS BY HARRY BATEMAN* 


. Question 14943. Educational Times (2) 1, 98-100 (1902). 
. Question 15119. Educational Times (2) 3, 110-111 (1903). 
. Question 15221. Educational Times (2) 4, 88 (1903). 
. The determination of curves satisfying given conditions. Proc. C. P. S. 12, 163-171 (1903). 
. Question 15440. Educational Times (2) 5, 68 (1904). 
. Question 15388. Educational Times (2) 5, 105-106 (1904). 
. The solution of partial differential equations by means of definite integrals. Proc. L. M.S. (2) 1,451-458 
(1904), 
8. Certain definite integrals and expansions connected with the Legendre and Bessel functions. Mess. (2) 
33, 182-188 (1904). 
9. A generalisation of the Legendre polynomial. Proc. L. M. S. (2) 3, 111-123 (1905). 
10. The Weddle quartic surface. Proc. L. M. S. (2) 3, 225-238 (1905). 
11. The correspondence of Brook Taylor. Bibliotheca Math. (3) 7, 367-371 (1906). 
12. Note on the solution of linear differential equations by means of definite integrals. Mess, (2) 35, 140— 
141 (1906). 
13. The theory of integral equations. Proc. L. M. S. (2) 4, 90-115 (1906). 
14. On the inversion of a definite integral. Proc. L. M. S. (2) 4, 461-498 (1906). 
15. Sur l’équation de Fredholm. Darb. Bull. (2) 30, 264-270 (1906). 
16. A class of integral equations. Trans. C. P. S. 20, 233-252 (1906). 
17. A type of hyperelliptic curve and the transformations connected with it. Quart. J. Math. 37, 277-286 
(1906). 
18. On an expansion of an arbitrary function in a series of Bessel functions. Mess. (2) 36, 31-37 (1906). 
19. On definite functions. Mess. (2) 37, 91-95 (1907). 


SDA 


* This bibliography was prepared by Joan Bateman. The following abbreviations are used: Bull. 
A. M. S.=Bulletin of the American Mathematical Society; Mess. = Messenger of Mathematics; Proc. 
C. P. S. = Proceedings of the Cambridge Philosophical Society ; Proc. L. M. S. = Proceedings of the London 
Mathematical Society; Proc. N. A. S.= Proceedings of the National Academy of Sciences; Trans. A. M. S. 
= Transactions of the American Mathematical Society; Trans. C. P. S.=Transactions of the Cambridge 
Philosophical Society. 
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20. The inversion of a definite integral. Math. Ann. 63, 525-548 (1907). 

21. The application of integral equations to the determination of expansions in series of oscillating functions. 
Trans. C. P. S. 20, 281-290 (1907). 

22. The reality of the roots of certain transcendental equations occurring in the theory of integral equations. 
Trans. C. P. S. 20, 371-382 (1907). 

23. (H. B. and D. M. Y. Sommerville). Question 16009. Educational Times (2) 11, 57-61 (1907). 

24. Some geometrical theorems occurring in hydrostatics. Mess. (2) 37, 119-123 (1907). 

25. A formula for the solving function of a certain integral equation of the second kind. Mess. (2) 37, 179- 
187 (1908). 

26. Notes on integral equations. 1. The integral equation of the first kind. Mess. (2) 38, 8-13 (1908). 

27. Notes on integral equations. 11. The method of least squares. Mess. (2) 38, 70-76 (1908). 

28. On the application of integral equations to the determination of upper and lower limits to the value of a 
double integral. Trans. C. P. S. 21, 123-128 (1908). 

29. On essentially positive double integrals and the part which they play in the theory of integral equations. 
British Ass. Rep. (Leicester) 77, 447-449 (1908). 

30. Question 16090. Educational Times (2) 13, 72-74, 91 (1908). 

31. The tangent plane which can be drawn to an algebraic surface from a multiple line. Arch. der Math. und 
Phys. (3) 13, 48-51 (1908). 

32. A method of calculating the number of degrees of freedom of a molecule among which the partition of 
energy is governed by the principal temperature. Manchester Mem. and Proc. 53, 1—9 (1908). 

33. The solution of linear differential equations by means of definite integrals. Trans. C. P. S. 21, 171-196 
(1909). 

34, Notes on integral equations. 111. The homogeneous integral equation of the first kind. Mess. (2) 39, e 
6-19 (1909). 

35. The conformal transformations of a space of four dimensions and their applications to geometric optics. E 
Proc. L. M. S. (2) 7, 70-89 (1909); British Ass. Rep. (Dublin) 78, 627-629 (1909). 

36. The reflexion of light at an ideal plane mirror moving with a uniform velocity of translation. Phil. Mag. 
(6) 18, 890-895 (1909). 

37. The solution of a system of differential equations occurring in the theory of radio-active transformations. : 
Proc. C. P. S. 15, 423-427 (1910). = 

38. The linear difference equation of the third order and a generalisation of a continued fraction. Quart. J. 4 
Math, 41, 302-308 (1910). 

39. The solution of the integral equation connecting the velocity of propagation uf an earthquake-wave in the 
interior of the earth with the times which the disturbance takes to travel to the different stations on the 3 
earth's surface. Phil. Mag. (6) 19, 576-587 (1910); Physik. Zs. 11, 96-99 (1910). = 

40. Notes on integral equations. 1V. The expansion theorems and the integral equation of the first kind. f 
Mess. (2) 39, 129-135 (1910). 

41. Notes on integral equations. V. Integral equations with variable limits. Mess. (2) 39, 173-178 (1910). 

42. Notes on integral equations. V1. The homogeneous integral equation of the first kind. Mess. (2) 39, 182- 
191 (1910). 

43. Report on the history and present state of the theory of integral equations. British Ass. Rep. (Sheffield) 
80, 345-424 (1910). 

44, The determination of solutions of the equation of wave motion involving an arbitrary function of three 
variables which satisfies a partial differential equation. Trans. C. P. S. 21, 257-280 (1910). 

45. Question 16215. Educational Times (2) 18, 86-87 (1910). 

46. A system of circles derived from a cubic space curve and the properties of a certain configuration of fifteen 
lines. Mess. (2) 40, 81-87 (1910). 

47. Kummer’s quartic surface as a wave surface. Proc. L. M. S. (2) 8, 375-382 (1910). 

48. The physical aspect of time. Manchester Soc. Mem. 54, 13 p. (1910). 

49. Correction to Mr. H. Bateman's paper on the Reflexion of light at an ideal plane mirror moving with a 
uniform velocity of translation. Phil. Mag. (6) 19, 824 (1910). 

50. Elementare elektronensysteme. (Elementary systems of electrons.) Physik. Zs. 11, 318-320 (1910). 

51. The transformation of the electrodynamical equations. Proc. L. M. S. (2) 8, 223-264 (1910). 

52. The transformation of coordinates which can be used to transform one physical problem into another. 
Proc. L. M, S, (2) 8, 469-488 (1910). 

53. The relation between electromagnetism and geometry. Phil. Mag. (6) 20, 623-628 (1910). 

54. On the probability distribution of a-particles. Phil. Mag. (6) 20, 704—707 (1910). 
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CONTRIBUTIONS TO THE PROBLEM OF APPROXIMATION 
OF EQUIDISTANT DATA BY ANALYTIC FUNCTIONS* 


PART B—ON THE PROBLEM OF OSCULATORY INTERPOLATION. 
A SECOND CLASS OF ANALYTIC APPROXIMATION FORMULAE 


BY 


I. J. SCHOENBERG 
University of Pennsylvania and Ballistic Research Laboratories, Aberdeen Proving Ground 


Introduction. The present second part of the paper has two objectives. Firstly, 
we wish to carry further the important actuarial work on the subject of osculatory 
interpolation (Chapters I and II). Secondly, we construct even analytic functions 
L(x), of extremely fast damping rate, such that the interpolation formula of cardinal 


type 


F(x) = yL(x — ») (1) 
reproduces polynomials of a certain degree and reduces to a smoothing formula for 
integral values of the variable x (Chapter III). This second problem is found to be 
intimately connected with the subject of osculatory interpolation. 

A preliminary remark concerning our notation is necessary. In Part A, Section 
2.21 we described various characteristic properties (or “type characteristics”) of a 
polynomial interpolation formula of the form (1), such as: (i) The degree m of the 
composite polynomial function L(x); (ii) its class C*, i.e., order of contact is yp; 
(iii) the highest degree k of polynomials for which the formula (1) is exact; (iv) the 
span s of the even polynomial function L(x). For convenience we propose to summa- 
rize all these statements by saying that (1) is a formula of type! 


* Received Oct. 18, 1945. Part A of this paper appeared in this Quarterly 4, 45-99 (1946). 

1 The connection of these types characteristics with the notation as used by Greville in his paper 
The general theory of osculatory inter polation, Trans. Actuar. Soc. Amer., 45, pp. 202-265 (1944), especially 
pp. 210-211, is as follows: The first three symbols D™, C*, E*, require no further comment since they are 
identical respectively with the characteristics 4, 1, and 6 of Greville’s classification, pp. 210-211. There 
remain three further characteristics to be discussed: (i) Whether the formula (1) is an “end-point” or 
“mid-point” formula. This point is of importance if (1) is written in terms of certra! differences, since it, 
then reduces to either the Everett or else the Steffensen form. The following statement is obvious: The 
formula (1) is an “end point” or “mid-point” formula depending on whether the span s is even or odd. 

(ii) Greville’s adjectives “ordinary” and “modified” agree respectively with our “ordinary” and 
“smoothing.” 

(iii) The highest order d of differences involved (explicitly or implicitly). We start with the following 
question: Let x be given. How many ordinates y, enter into the computation of F(x) by (1)? Assuming that 
L(x) is continuous, hence =0, at the end point x=s/2 of its span, we have L(x—n) 40, as long as n in 
such that 
|x —n| <s/2. 


This inequality is found to be equivalent to 
: <an< (*) 
——+ —+x. 


Let s =2¢ be even (end-point formula) and let now x be anywhere within 0Sx 31. By (*) F(x) then re- 
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D™, C4, E*, s. (2) 


As an instance we may describe the k-point central interpolation formula of Part A, 
Section 2.121 as a formula of type? 


pe .. if k is even 


E*|, s= k. 
C—' if k is odd 


In Chapter II we construct a class of ordinary interpolation formulae and two 
classes of smoothing interpolation formulae. These classes by no means describe all 
possible osculatory interpolation formulae. Furthermore, a number of interesting 
problems concerning remainder terms and orders of approximations await solution. 
No attempt has been made to see which of the numerous formulae tabulated by 
Greville, loc. cit., are contained in the three classes of formulae developed in Chapter 
II. The essential progress made in this direction may perhaps be briefly described as 
follows. The construction of an interpolation formula usually requires the solution of 
a more or less complicated system of linear equations, unless, as in Lagrange’s for- 
mula, the basic interpolating functions are obvious from the start. These systems of 
equations are especially troublesome if one wishes to construct an osculatory inter- 
polation formula of any general class. As Greville correctly points out, loc. cit., pp. 
255-256, the mere agreement between the number of unknowns with the number of 
equations which they should satisfy, will, by itself, never prove the existence of a 
solution. Basically, our parametric representation of spline curves of order k (Part A, 
Section 3.15, Theorem 5) circumvents this difficulty. 

An example which illustrates the operation of this representation is as follows. 
Let F(x) be defined as equal to 0 for x $0, as well as for x24. We propose to complete 
the definition of F(x) in the range 0S x4 by four cubic arcs joining at x =1, 2, 3, in 
such a way that F(x) be of class C”’ for all real x. Of course, we are not interested in 
the obvious but trivial solution F(x) =0. Let us now count the available parameters 
and the number of conditions. The 4 cubic arcs furnish 4-4=16 parameters. The sec- 
ond order contact requirements at x =0, 1, 2, 3, 4 lead to a system of 3-5=15 homo- 
geneous equations. The solution of a homogeneous system of 15 equations in 16 
unknowns depends on anything from 1 to 16 arbitrary parameters, depending on the 
rank bf the system. Our Theorem 5, for k=4, furnishes immediately the one-parame- 
ter solution 
F(x) = ¢-M,(x — 2) (3) 
the graph of which is given in Part A, Section 3.13. Again Theorem 5 will easily show 
that this is the most general solution of the problem. We see how this complicated 
system of 15 equations in 16 unknowns is explicitly solved by (3). As a variation of 
the problem, let us now define F(x) to be equal to 0 for x £0, as well as for x23, and 
let us propose now to bridge this gap by 3 cubic arcs giving a F(x) of class C’. Now we 
find that the problem amounts to a system of 12 homogeneous equations in 12 un- 
knowns. This tells us precisely nothing. Again by Theorem 5 we can readily show that 


quires all y, such that —7o<n<o-+1 that is s=2s consecutive ordinates. Let s=20+1 be odd (mid-point 
formula) and let x be anywhere within —}Sx}. Again by (*) F(x) now requires all y, such that 
—o—1<n<o-+1, hence again s=20¢+1 consecutive ordinates. We have therefore proved the following: 
The highest order d of differences involved is clways related with the span s by the relation s=d+1. 

2 The symbol C~ is to indicate the class of piecewise continuous functions. 
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the trivial solution F(x) =0 is the only solution. These considerations generalize and 
allow to characterize our basic functions 
1 k k-1 
M(x) = ———_i x, , 
(k — 1)! 
up to a multiplicative constant and a shift along the x-axis, as follows:* Let F(x) 
be =0 for x $0, as well as for x2=n, where n 1s a positive integer. We wish to complete 
the definition of F(x) by a succession of n arcs, of degree k—1, joining atx=1,2,---, 
n—1 such as to furnish a F(x) of C*-?. Then n=k is the smallest value of n for which 
this can be done in a non-trivial way and for this minimal value n =k the gap is bridged by 


F(x) = ¢-M,(x — k/2) 


and in no other way. 
The reader who is mainly interested in the numerical applications may pass di- 
rectly from here to the Appendix where the use of the tables is fully explained and one 


example is worked out. 


I. THE COSINE POLYNOMIALS :(u) AND CERTAIN RELATED SETS OF POLYNOMIALS 


In the present chapter we propose to study further properties of the cosine poly- 
nomials 
¢x(u) = >> M,(n) cos nu (1) 
which were mentioned in Part A, sections 3.14 and 4.1 (for t=0). By Part A, section 
4.1, formula (6) (for t=0) we may also write 


o.(u) = vi(u + 


and therefore ssi 
ox.(u) = (2 sin u/2)* 
+ 
1.1. Expression of ¢,(u) in terms of rational polynomials. We introduce two new 


sets of periodic functions defined by 


(2) 


= 1 
px(u) = (2 sin u/2)*- ae (3) 
ox(u) = (2 sin u/2)*- . (3’) 
+ 
A comparison with (2) shows that : 
px(u) if k is even, 
if is odd. (4) 


3 Problems of this kind concerning polygonal lines of a certain degree and class are of importance for 
the theory of formulae of mechanical quadratures. The author expects to discuss this connection else- 


where. 
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By differentiation of (3) and also (3’) we readily find the recurrence relations 


u 
= Cos > px(u) — (u), (5) 
u 
= ry sin OK (u). (5’) 


These may be used in a progressive computation of p, and o, if we start with 


p2o(u) = 1, o,(u) = 1. (6) 
We prefer, however, to express p;(u) and o,(u) as polynomials in the variable 
x = cos (u/2) (7) 
by means of 
pi(u) = Ux_2(cos u/2), ox(u) = Vi_s(cos u/2). (8) 


Substituting into (5), and (5’) respectively we find that the two sequences of polyno- 
mials U(x), V.(x), both of exact degree n, satisfy the recurrence relations 


1 
Visi(x) = + (1 — (x), (9) 
1 
Viri(x) = xVi(x) + (1 — (x), (9%) 
with initial values which by (6) and (8) are 
U(x) = 1, Vo(x) = 1. (10) 
A simple calculation now gives 
U(x) = x, Vi(x) = x, ) 
1 1 
U(x) = (1 + = (1+ 2’), 
= + 29 = — (Sx + 29) 
x) = —(2x+ x), =—(5*+x 
1 
= — (2 + 11x? + 2x4), Va(x) = — + + x4) 
15 24 
1 1 
U;(x) = (17x + + 225), V;(x) = (61x + + x5), 


We record as a lemma the following properties which are readily established by in- 
duction. 


LemMA 1. U;(x) and V;(x) are polynomials of exact degree k which are even or odd 
according as k is even or odd. The coefficients of their highest terms are positive. Also 


U1) = Vi(1) = 1, 1) = 1) = (— (12) 
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In view of (4) and (8) we find the following expression of $;,(u) in terms of our new 


polynomials: 
U;-2(x) if kis even, 


if kis odd. 


o.(u) = { (13) 


This shows that the even polynomials U2,, V2, are of special interest. 
1.2. The zeros of the polynomials U; and V;. We propose to prove the following 


proposition: 
LEMMA 2. The zeros of the even polynomials U2,(x) and V2,(x) are all simple and 
purely imaginary. 


We carry through the proof for U2,(x) only since the proof for V2,(x) is entirely 
similar. In order to deal with real zeros, we define a new sequence of polynomials 


ux(x) by 
= (k = 0,1,---). (14) 


These new polynomials are also real and satisfy a recurrence relation which in view 
of (9) is readily found to be 


= xu,(x) — (1 + x*)u;! (x). (15) 


Froin (11) we find 
u(x) = 1, u;(x) = x, uo(x) = 3(2x? — 1), us(x) = 3(x? — 2x),---. 
In view of (14) it obviously suffices to show that the zeros of u,(x) ar real and simple, 


while those of u2,(x) are also different from zero. This is readily done by induction 
as follows. Let k =2v be even and Jet us assume that the k zeros of u,(x) are 


and therefore simple. This, and the fact that u,(x) has a highest term of positive co- 
efficient (Lemma 1 and (14)), imply that 


ug (é,) > 0, 


and that the sequence of values of u/ (x), at the k roots (16), alternate in sign. By 
(15) we therefore find 


<0 


and that the values of ux.4:(x),at the k roots (16), alternate in sign. Since u44:(0) =0, 
we conclude that u44:(x) has v positive and v negative zeros which must therefore be 


simple. 
Let now k=2v+1 be odd and let u have the simple zeros 


— ~ <&). (17) 


Now we conclude as before that u,,,(£,) <0 and that the values of u;4;(x), at the k 
roots (17), alternate in sign. Again the conclusion is that u44:(x) has simple real roots 
none of which vanishes. This proves the theorem by complete induction. 

1.3. A few corollaries. In this last section of the present chapter we prove several 
auxiliary propositions which will be used in the next chapter in the derivation of inter- 
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polation formulae of various kinds. These propositions represent the solutions of the 
algebraic problems arising by the Fourier integral transformation of the problems of 
the construction of those interpolation formulae. 


LemMaA 3. Let k=2v be even. We can determine uniquely an even polynomial P;,(x), 
of degree k, and an odd polynomial P;,_,(x), of degree k—1, satisfying the identity 


+ = 1. (18) 


Likewise polynomials Q;(x) and Qy-1(x), even and odd respectively, exist uniquely such 
as to satisfy 


+ = 2. (19) 


We wish to show first that U; and Uj4; have no common zeros. Indeed a common 
zero x of U;, and U;41 would, by (9), be a zero of 


(1 — (x). 


Since by (12) «# +1, x must be a zero of U;/ (x). But this contradicts our Lemma 2 
to the effect that U;.(x) has only simple zeros. The polynomials U;(x), Ux4:(x) having 
no common divisors, the identity (18) is assured by the elementary theory of the great- 
est common divisor of two polynomials. We now show that P;(x) is even and P;_;(x) 
is odd as follows. Replacing x by —x in (18) we find 


Ui(x)Pi(— x) — Unqi(x)Pia(— x) = 1. 
Since our polynomials P;, P;,_; are uniquely defined by (18) we find 
= x), = — Pia(— 2), 


which prove our statement. An identical reasoning proves the existence of the poly- 
nomials Q; and Q;41 satisfying (19). 

The polynomials P; and P,4; are easily determined for low values of k by the 
method of indeterminate coefficients. Thus for k=2 by (11), 


U2(x) = (1 + 22*)/3, Us(x) = (2x + x*)/3, 
from which we find 
P.(x) = 2x? + 3, P(x) = — 4x, 

satisfying the identity 

+ U3(x)Pi(x) = 1. (20) 
Likewise for k=4 we have by (11) 

Va(x) = (5 + 18x? + x4)/24, Ve(x) = (61x + 58x* + 2°)/120. 
The corresponding polynomials Q,, Qs are found to be 
Qu(x) = (3648 + 47892? + 83x4)/760,  Qs(x) = — (1469x + 83x3)/152. 

They satisfy the identity 

Va(x)Qa(x) + Vs(x)Qs(x) = 1. (21) 
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The identities (18), (19) will later be used in the following form. Again for an 
even k, but replacing k by k—2, we get by (18) and (8) 
+ = 1, (keven, x = cos u/2). 
Likewise for an odd k, but replacing k by k—1, we obtain from (19) and (8) the 


identity 
ox(U)Qn-1(4) + on41(#)Oxn-2(k) = 1, (k odd, x = cos u/2). 


The even polynomials x~'Px-3, and Qy-1, may now be expressed in 


powers of 
1 — x? = (sin 


We have therefore proved the following: 
LemMA 4. We can find constants a,, a; , b,, b/ such as to satisfy the following two 
identities : 
For an even k 
px(u) {ao — sin u/2)? + a4(2 sin u/2)* — -- + + ayo(2 sin u/2)*-*} 
+ prsi(u) {ad — ag (2 sin u/2)? + af (2 sin u/2)4 — --- 
ag_4(2 sin u/2)*-4} (2 cos u/2) 


1, (22) 
and for k odd 


ox(u)\bo — sin u/2)? + b4(2 sin u/2)* — + sin u/2)*—} 
+ onsi(u) {bg — bs (2 sin u/2)? + bf (2 sin u/2)4 — --- 
F sin u/2)*-#}(2 cos u/2) = 1. (22’) 


As examples we mention that the identities (20) and (21) become on passing to the 
variable u 


pa(u){5 — 3(2 sin u/2)*} + ps(u){— 2}(2 cos u/2) = 1 (23) 


and 


(2 sin + —— (2 sin w/2 ‘ 
760 12160 7 


19 
sin =) cos = 1. (24) 


The last proposition which we wish to derive here concerns the expansion of 1/¢:(u) 
in ascending powers of the variable 


s = sin? u/2 = 1 — cos? u/2 = 1 — x? (25) 

Let us assume for the moment that & is even. Then by (13) 
= Ux-2(x), (k even). (26) 
Now U;-2(x) is an even polynomial which, by Lemma 2, has purely imaginary zeros. 
Being an even polynomial, U;-2(x) may be expressed as a polynomial U*(s) in the 


variable 
s=1— 2? (25’) 
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of degree x= (k—2)/2. This change of variable transforms the purely imaginary zeros 
of Ux_2(x) into the zeros 

2, °° * Oy (x = (k — 2)/2) 


of U*(s) which, by (25’), must all be positive and greater than 1. Finally, since 
U,-2(1) = U*(0) =1, we have the identity 


An entirely similar identity is derived for an odd k by repeating our arguments for 


ox(u) = Ve_s(2), 
instead of (26). 
This establishes the following 


Lemma 5. The reciprocal of the cosine polynomial $;(u) admits of an expansion 
1 “ed n 
= (2 sin u/2)" (28) 
n=0 


which converges for all real values of u and where the coefficients are positive rational 
numbers 


(k) 


>¢ = 0,1, 2,---). (29) 

Indeed, in view of (27), the expansion (28) may be obtained as 
1 
= 
I] ( ay a 


which reduces to (28), in view of (25). In conclusion we notice the following conse- 
quences of the identity (28). Since 


+++ ) con (48)" 


n=0 


s 
2 


2 sin u/2\* 
vi(u) = 


we have in the neighborhood of the origin u=0 


where 


2 si 2\* 
= + u*-(regular function of (30) 
u 
On multiplying (28) by ¢4(u) we therefore have 


2sinu/2\* <= " i 
( ) Con (2 sin = 1+ (regular function) 
u n=0 


and also 


sin «t+ -(regular function). (31) 


OS2n<k 


& 
(2 7) 
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It is of special interest to point out that if 


4 
Sk,m(u) = : con (2 sin u/2). 
u 


n=0 
then 
1+ u?”-(regular function) if 2m< k 
1+ u*-(regularfunction) if 2m—2< k S 2m. 
As an illustration we find for k=6 by (13), and (11), and (25) 
1 1 15 30 1 
be(u)  Us(x) 2+ 1122+ 2x4 30 — 30s + 4s? 2 
15 
whence 
1 is 
15 
or 
= 1+ — (2 sin u/2)* + (2 sin u/2)* + (33) 
= — (2 sin u/2)? + —(2sin u HHA, 
4 240 
The relations (32) now become (for k=6, m=1, 2, 3) 
2 sin u/2\° 
= 1+ u*-(regular function), 
u 


Il 


1 + (regular function), 


{1 + (2 sin w/a 


(- sin {i 1 (2 sin u/2)? + 13 (2 si 
a — n u/2)? + — (2 sin 
sin 4 in u 


1 + u®- (regular function). 


In the next chapter we shall need the numerical values of the coefficients 
(k) 
Go for 2n< k. (34) 
It is of interest then to point out that these coefficients (34) may also be otherwise 
computed as coefficients of a simple generating function. Indeed, from (30) it is clear 
that the coefficients (34) will not change if on the left-hand side of (28) we replace 
¢,(u) by the first term on the right-hand side of (30). That means, if 


> dy, (2 sin u/2)” (35) 


2 sin u/2 por 


then 
(k) (k) 


Can, = day (2n < k). (36) 
However, the coefficients of (35) are readily determined. Indeed, if we set 
vy=2sinu/2 or u = 2 arcsin 0/2, (37) 


then (35) becomes 
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2arcsiny/2\* <= 
= (38) 
v n=0 
Since 
2 arcsin 0/2 3 ¥ 1 
v 23 4 2-45 16 2-4-6 7 64 


we find the expansion (38) by raising (39) to the k-th power. Thus 


k 5k? + 22k 
Cat 
24 5760 


Since all coefficients of (39) are positive, it is clear that the coefficients of (38) are 
likewise positive. This however does not imply the positivity of the coefficients of 
(28) beyond the kth term. For k =6 the values (46) agree with the coefficients of (33). 


(40) 


II. POLYNOMIAL INTERPOLATION FORMULAE 


In this chapter we wish to apply our general Theorem 2 of Part A, section 2.23 
and our Lemmas 4 and 5 of the last section in deriving three distinct classes of poly- 
nomial interpolation formulae for each value of the positive integer k. The formulae 
of the first class (Theorem 1 below) are of the ordinary kind (see Part A, section 2.21 a 


and b), and of the type 


if kis even 
Ce*, s= 


2k—1 if kis odd. 


The existence of ordinary interpolation formulae of degree k and class k—2 was pre- 
viously conjectured by Mr. Greville who verified their existence up to and including 
k=6. (See Greville, loc. cit., pp. 212-213.) The formulae of the second class are 
smoothing interpolation formulae (Theorem 2 below). For a given integral k and 
each integral m, such that 0 $2m—2<zk, a formula is derived which is of the type 


D*", Cc, Emin(2m—1,k-1) | s= k 2m 2. 
These formulae are derived from an ordinary interpolation formula of type 


E*, $= ©, 
discussed in Part A. 

The formulae of the third and last class are again smoothing interpolation formulae 
(Theorem 3 below). While in the second class the degree D*~! and the “order of con- 
tact” C*~? were fixed, while the degree of exactness E*"~! and the span s=k+2m—2 
increased apace, in the present class the span s=k is constant. More precisely, a 
formula is derived for each m such that 0<2m<k—1 which is of type 


Ch-2m, E21, s= b. 
These formulae are derived from the formula of ordinary k-point central interpolation 
in a manner somewhat reminiscent of Mr. Jenkins’ original procedure. 
2.1. Ordinary polynomial interpolation formulae of the Jenkins-Greville type. 
We are returning to our basic functions 


1 /2 sin u/2\* 
M(x) = e“*du (1) 


| 
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and wish to show that the osculatory interpolation formulae of the type investigated 
by Jenkins and Greville may be readily derived in terms of these functions. We shall 


use the operational symbol ¢ to mean 


of(x) = f(x + 3) + f(x — 2). 


THEOREM 1. We define the basic polynomial function L(x) by the following two for- 


mulae according to the parity of k: 
L(x) = aoM,(x) + + +++ + 

+ ag oM + ag + + (R even) 
L(x) = boMi( x) + + + (x) 

where the numerical constants a,, a; , b,, b/ are those defined in Lemma 4. Then 


F(x) = ynL(x — n) 


is an ordinary polynomial interpolation formula of type 
—2 if kis even, 
D', Ce, EM, s= 
2k—1 if kis odd. 
Indeed, we notice that 
= e~iu/2) giuz = 2isin u/2e“*, 


= (eiu/? + tul2) gius = 2 cos u/2e*“*, 


Let k now be even, hence L(x) defined by (2), and let 


L(x) = — 
We evidently obtain this integral representation by performing the operation 
+ + + ay—25*-? 
on the relation (1) and add to it the result of performing the operation 
ago + +--+ + 
on (1), with k replaced by &+1. In view of (6) we have 
= (27 sin u/2)’e'"* 
= (2i sin u/2)"(2 cos u/2)e™?, 
and therefore 


2 sin u/2\* . 
g(u) = {ao — ao(2 sin u/2)? + +++ + sin u/2)*-*} 
u 


2 sin u/2 
(a 


k+1 
) {aj — af(2sin u/2)?+--- 


F aj_4(2 sin u/2)*-*} (2 cos u/2). 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 
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We now turn to Theorem 2 of Part A, section 2.23, which states that (4) is an ordinary 
interpolation formula if and only if the following identity holds: 


g(ut+ 2m) = 1. (8) 
It should be noticed now that both expressions in (7) contained within braces are 
periodic functions of period 27. Since k is even we find that }-g(n+2zv) is identical 
with the left-hand side of our relation I (22). This proves (8) for even k. A precisely 
similar reasoning for odd & will show that }>g(u+2zv) is identical with the left-hand 
side of I (22’). 
There remains the problem of showing that (4) is of the type as stated in the Theo- 
rem. Since L(x) is by (2), (3), a linear combination of functions of the form 


M(x +n), + 3+), 


it is clear that L(x) is a polynomial line of degree k, of class C*-*, with discontinuities 
at x=n, or x=n+1/2, according to whether & is even or odd. Finally (4) is exact for 
the degree k—1, again by Theorem 2 of Part A, section 2.23. There remains the 
discussion of the span s of L(x). Now the span of M;(x) is =k (see Part A, section 
3.13) and therefore the span of 6”M;(x) is equal to k+v, while the span of 06’M;(x) 
is equal to k+v+1. Now it is immediately verified that the two terms of (2) involving 
5*-? and 76*~‘ are both of span 2k —2, while the similar two terms of (3) are both of 
span 2k—1. This completes a proof of the Theorem. 

As illustrations we mention that the identities I(23) and I(24) corresponding to 
the cases k=4 and k=5 give rise to the basic functions 


L(x) = 5M4(x) + 369M,(x) — 20M;(x) (9) 
and 
L(x) = + s(x) + 54M 
19 760 12160 
194 
oM,(x) — (10) 


These two basic functions give rise to ordinary interpolation formulae (4) which are 
of the types 
D,C?, E,s=6 and C*, Et, s =9, (11) 


respectively. 
Incidentally, the characteristic function of (9) is, by I(23), 


g(u) = (s — 2sin? —4 cs u/2 
u 


2 sin u/2\4 sin “ 
g(u) = (4 +cosu—4 =*). 
u 


u 


or 


This agrees with our formula (11’’) of Part A, section 2.122, as in fact the basic func- 
tion (9) is identical with Jenkins’ function there described by formula (11). 


ay 
| 
q 
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In concluding we wish to mention the numerical results for k =6. In this case we 


need the identity 
Us(x)Pa(x) + Us(x)Ps(x) = 1. (12) 


By (11) we have 


1 1 
= (2+ 2x4),  Us(x) = (17x + + 2x5). 


By indeterminate coefficients we find 


1§ 115 27 285 27 
Px) = —-— — — 2, 


from which, on passing to the variable I(25), the identity (12) becomes 


14-28 448 


339 
+ p7(u) {- a + = (2 sin w/2y\ (2 cos u/2) = 1. 


The basic function corresponding to k=6 is therefore 


giving rise to an ordinary interpolation formula of type 
D*, E® and s = 10. 


2.2. A first class of smoothing interpolation formulae derived from an ordinary 
interpolation formula of type D*—', C*-*. We start by recalling an ordinary polyno- 
mial interpolation formula derived in Part A, section 4.2. Indeed, the formula (9) of 
that section furnishes, for t=0, the following polynomial basic function 


1 
Lz) = — (14) 
px (u) 
The corresponding formula 
F(x) = D> yali(x — n) (15) 


is, as we know, an ordinary polynomial interpolation formula of type 
{* if 3, 
k if k=1, 2. 


We turn now to the expansion I(28) of Lemma 5, substituting I(28) into the integral 
(14) and integrating term-wise we obtain the expansion 


Li(x) = — 6s 8 +o, 8 Mix) (kz 3). (17) 


: 
n=—0o 
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On comparing the present interpolation formula (15) with the formula (4) of the 
Jenkins-Greville type, we notice, by (5) and (16), that they are both of class C*-? 
and that they are both exact for the degree k —1. The degree of (15) is lower by 1 than 
the degree of (4). This reduction of the degree to the lowest possible value k—1, for 
a formula of class C*-*, was achieved at the price of having an infinite span. The 
infinite span of (15) clearly disqualifies this interpolation formula as far as numerical 
purposes are concerned. 

We now turn to the partial sums of the series (17). They will yield smoothing in- 
terpolation formulae of considerable practical importance. Indeed, let 

m—1 (k) 


= Mi(x) — Mi(x) + (— 1)” Mi(x), (2m — 2 < k). (18) 
The characteristic function of this basic function is identical with the left-hand side 
of the identity 1(32). In view of our Theorem 2 of Part A, this identity 1(32) proves 
that (18) is the basic function of a smoothing interpolation formula which is exact 
for the degree equal to min(2m—1, k—1). It is, moreover, visibly of degree k—1, of 
class C*~*, and of span s=k-+-2m—2. One further important point is in need of proof, 
namely that the formulae based on (18) actually do smooth any given sequence (see 
Definition b of Part A, section 2.2). This will readily follow from Lemma 5. Indeed 
the characteristic function @m,,(u) of the formula 


F(n) = — (19) 
is, by Theorem 2, Part A, given by 


= > £k,m(u + 


By 1I(32) and I(28) we now have 
= > + 2xv) = 6s. (2 sin 
< (2 sin u/2)"=1, (0<u<2n). (20) 


n=0 


Since obviously @m,.(u) >0, for all u, we see that (19) is indeed a smoothing formula 
according to our definition. Recalling the relations I(36), 1(38), we may therefore 
state the following Theorem: 


THEOREM 2. Let k be a positive integer and m an integer such that0<2m<k+2. Let 
the positive rational numbers dS be defined by the expansion 


v n=0 

Then 

= Milx) — dy 8 + dy 8 Mila) — + (— 1)” (22) 


gives rise to a smoothing interpolation formula 


q <? ed 
‘ae 
ig 
| 
a 
= 
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F(x) = ynLe.m(x — m) (23) 
of type 
D*-|, Emin(2m—1,k—-1) s= k 2m (24) 


Moreover, the formula (23) preserves the degree k—1. (See Part A, section 2.21, Defini- 
tion d.) 

If k is fixed and m increases then the smoothing power of our formula (23) decreases 
according to our definition of Part A, section 1.12, Definition 2. 


The last statement concerning the decreasing smoothing power of (23) follows 
from (20), since ¢k,(“) increases strictly as m increases while u is constant 


(0<u<2z). 
Notice, by (24), how on increasing m by one unit both the degree of exactness, as 


well as the span, increase by two units. 
As illustration we find from I(40) that 


k 
Li2(x) = Mi(x) — 5°Mi(x), (k2 4), (25) 
yields a smoothing formula of type 
p>, E*, s=k+2. (26) 


The characteristic function of (25) is 


(i =) (27) 
= ( ) sin 


2 sin u/2\* k k 
= {1 + — — — cos 
u 12 12 


For k=4 this function g42(u) agrees with the integrand of our formula (12’) of 
Part A, section 2.123. Also Mr. Jenkins’ basic function, as given by formula (12) of 
Part A, section 2.123, may be derived by working out the various polynomial expres- 


or 


sions of 
L4.(x) = 4( x) (28) 


from the explicit expressions of M,(x) (see Part A, section 3.13, (14)). 
Likewise, by 1(40) 


k + 22) 
= M,(x) — — + 54M (k 2 6), (29) 
24 5760 
yields a smoothing formula of type 


2.3. A second class of smoothing interpolation formulae derived from the ordi- 
nary k-point central interpolation formula. Among the smoothing interpolation for- 
mulae (23) described by Theorem 2 the one of most interest is obtained by letting m 


| 
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assume its largest value. If k is even, m is maximal if 2n—2=k—2 or m=k/2. If k 
is odd, m is maximal if 2m—2=k—1 or m=(k+1)/2. In either case 


[=]. (31) 


where [x] represents the largest integer not exceeding x. The corresponding basic 
function (22) is 


x k 
=} {1 + ds°(2 sin + 


+ sin u/2)”\e (32) 


We recall that the smoothing interpolation formula based on this function is by (24) 
of the type 
ons, s=k+2u—2. (33) 


Indeed, the formula is exact for the degree k—1 because of 


2 sin u/2\* | 
= 1+ u*-(regular function). (34) 


An interesting counterpart to (32) is obtained as follows. An identity of the type 
(34) may also be obtained if in the expression within braces we replace 2 sin u/2 by u. 
Indeed, rational constants y{ may be determined such that 


2 sin u/2\* (k) 2 (k) 4 (k) 2 


= 1+ u*-(regular function). (35) 
LeMMA 6. The basic function 


Ty u(x) (= ) {1 + "te du (36) 


is identical, for all real values of x, with the basic function C,(x) of the k-point central 
interpolation method (see Part A, section 2.121). 


Notice first that by differentiation of 


1 sin u/2\* 
J u 


My (x) = (— 1)" =f < k). (37) 


we obtain 


Therefore the integral (36) may also be written as‘ 
4 Assuming (41) already established, we see by (38) and the relations 
= (S»Sk-1), (*) 


(see Part A, section 3.15, formula (23)) that we may express C;(x) as follows 


- 
| 
4 
= 
re 
‘ 
; 
ca. 
i 
\ 
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= k 
= Milt) — ME! (x) Me (2) — + 1)” (2). (38) 


Now M;(x) and all its derivatives are functions of span s=k. Therefore also 
T',,,(x) has the span s=k. Furthermore by (35), and Theorem 2 of Part A, we con- 
clude that 


is an interpolation formula of the following type characteristics: 
s=k, (40) 
These last two properties (40) allow us to show readily that 
= C,(x) for all real x. (41) 


Indeed, let k be even, k=2x. Let Po(x) be the polynomial of degree k—1 defined by 
the following k conditions 


P(—« +1) = P(—« +2) =--- =Po(—1)=0, = 1, 
Since (39) is exact for the degree k —1 we have the identity 
Po(x) = >> Po(n)Ti,(x — mn), for all real x. (43) 


We now restrict x to the range 
1. (44) 


Then we may write (43) as 
= D> Po(n)Ti,(x — n) 


n=—«+1 
since T;,,(x—m) =0 if |x—n| =x. In view of (42) this identity reduces to the single 
term, for n=0: 
Pox) =Tr,(x), 1), 
and therefore (41) holds for the range (44). Likewise, applying the formula (39) to 
the polynomial P(x), of degree k, defined by 


+2) =--- = P;(— 1) = 0, P,(0) = 1, P\(1) =--- = + 1) = 90, 


we find that (41) holds in the range 1 Sx 2, and so forth. Similar arguments obvi- 


ously apply, with obvious modifications, to the case of an odd k. 


The coefficients y¥ are the expansion coefficients of 


& 
(; sin =) van (45) 


Ce(x) = Malz) — + ala) (= Mele). 
This formula reveals at a glance the following fact: If k is even, then C,®”(x) (v=0, 1,2, +++ ) arecontinu- 
ous. If k is odd then C,’+)(x) (v=0, 1, 2, + + + ) are continuous. The author learned this property from 
Kingsland Camp, Notes on Interpolation, Trans. Actuar. Soc. Amer., 38, p. 22 (1937). 
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In N. E. Nérlund’s Differenzenrechnung, page 143, we find the expansion 


( 


sin ¢ (2r)! 
The coefficient D® is a polynomial in k of degree v. Nérlund’s Table 6 on page 460, 
loc. cit., lists these polynomials for v=0, 1, - - - , 6. We therefore have 
(k) 
————} = — 1)’ . 46 
€ sin =) (2v)! 2% (46) 
whence 
(k) 
(k) Dy 
47 
( ) (2v) (47) 
The first few values are 
k k(Sk + 2) 
(k) (k) (k) 
= 1, =— = ——_—_—_-- 48 
Yo v2 24" v4 5760 (48) 


The expansion coefficients of ¢/sin ¢ are positive (see Nérlund, loc. cit., Chapter II, 
sections 2, 3). Therefore the coefficients y are all positive. We shall use this fact later. 
In view of the results of our last section it seems natural to consider the partial 


sums 
m—1 (k) (2m—2) 


= Mix) — + +(— 1)” (2), 


of the sum (38). The properties of the interpolation formulae based on these functions 
are described by the following theorem: 


THEOREM 3. The formula 


k+1 


is a smoothing interpolation formula of the type 
Ck-2m, kh (51) 


The smoothing power of (50) decreases, as m increases from m=1 tom=p—1, until for 
m =p (50) reduces to the (ordinary) k-point central interpolation formula. 


Indeed, the characteristic function of (49) is 
2 sin u/2\* ~~ 
g(u) = ( ) (S2) 


From (45) we conclude that 


g(u) = 1 + (regular function) 


a 
iy 
| 
} 
a, 
= 
i 
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and therefore (50) is exact for the degree 2m—1, by Theorem 2 of Part A. The re- 
maining three characteristics 


D s= k, 
are evident on inspection of (49). There remains the investigation of the characteristic 
function of the corresponding smoothing formula 


F(n) = >> — v). (53) 
By Theorem 2, Part A, this characteristic function is 


xm(u) = >> g(u + 2nv). 


From (52) and I(3), we obtain 
= bu(u) + (2 sin u/2) +--+ + Sin 4/2)” (54) 


This expression is obviously positive for all values of u. For m=, however, we ob- 
tain the identity 

Xu(u) = 1 (55) 
since, by Lemma 6, we have before us an ordinary interpolation formula. Since (54) 
is a partial sum of the left-hand side of (55) we have therefore proved the inequalities 


0 < xm(u) < 1 (56) 


Then (53) is indeed a smoothing formula. The final statement of the Theorem is evi- 
dent from (54), since xm(u) increases with m. 

As illustrations we mention the following four special cases, two from each end 
of the range of values of m. 

(i) m=2. The formula based on 


k 

Tx,2(x) = M 24 (x) (k 4) (57) 

has the type 
E?, s = k. (57’) 

(ii) m=3. The formula based on 
k,3 & 24 k 5760 k ’ = ’ 

has the type 

E5, s= k. (58’) 


The values (48) were used. 
(iii) m=u—1. The formula based 


5 The formula (59) is especially instructive because we can observe very clearly how the addition 
to C;(x) of the extra term removes the crudest discontinuities of C,(x). Indeed, by the formula (*) of 
our preceding footnote we may write (59) as 
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a(x) = Cr(x) + (— (RS 4) 


has the type 
" if kis even 


C! if kis odd 


(iv) m=yu—2. The formula based on 


(k) (2u—2 (k) (2u—4) 


Te = Ce(x) + (— 1) (x) +(— 1) (x) (R26) (60 
has the type 
a if kis even 
D*-|, 


E*-5, s = k. (60’) 
C? if kis odd 


These formulae show clearly how an increase in “order of contact” is compensated 
by a corresponding loss in “reproductive power” and vice versa. 


III. A SECOND CLASS OF ANALYTIC INTERPOLATION FORMULAE 


In Part A, section 4.2, we described a class of ordinary analytic interpolation for- 
mulae of basic function 
1 t) 


t) = 


k=1,2,--+;#>0). 1 


These interpolation formulae are exact for the degree k —1. The basic function (1) as 
well as those of the smoothing interpolation formulae derived from it in Part A, sec- 
tion 4.3, dampen out like a descending exponential function. In the present last chap- 
ter we wish to construct smoothing analytic interpolation formulae of basic functions 
dampening out like 
exp (— c?x?), 

hence much more rapidly. In view of the development of section 2.2 it would seem 
fairly obvious how such formulae may be derived. We clearly need an analogue of 
Lemma 5 which we state asaconjecture: The reciprocal of $x(u, t) admits of an expan- 
sion 


2y-2 


= Ce(x) + (— 1) 
Let k be odd, hence 2u—2 =k—1 and therefore 
a(z) = Ce(x) + (— (59””) 


As seen from the graph of M4,(x), the corrective term is a step-function with discontinuities at x =n-+1/2 
whose values are proportional with the binomial coefficients of order k—1: (*F!). Their addition to Cy(x) 
offsets the discontinuities of C;(x) and turn it into a function (59’’) of class C!. If k is even, hence 2u—2 
=k—2, we have 

= + (— (59’””) 


As seen from the graph of M2(x), the corrective term is now an ordinary polygonal line with vertices at 
x=n, whose ordinates (at these vertices) are proportional to the binomial coefficients of order k—2: 
(*>*). Again, the superposition of this polygonal line on C;(x) offsets the corners of C;(x) and turns it 
into a function (59’’’) of class C?. The formulae (59’’’), (59’’) are especially convenient for constructing 
tables of these functions from existing tables of C;(x), i.e., tables of Lagrange interpolation coefficients. 
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sin 2/2)", (2) 


which converges for all real values of u and where the coefficients are all positive 


ca (t)>0, (n=0,1,2,---). (3) 


A proof of this conjecture would require a closer function-theoretic study of the 
entire periodic function ¢;(u, ¢) which has not been carried through as yet. 


Since 


t) = + t) 


(see Part A, section 4.1, formula (6)) we have 
t) = 4) + u*- (regular function). 


Therefore the expansion (2) agrees in its terms of order less than k with the similar 
terms of the expansion 


1 u . 


t) sin u/2 
Hence 
(k) (k) 
Can (t) = don (2), (0 s 2n < k). (5) 


The expansion (4) is readily determined and its coefficients are found to be positive 
as follows. We turn back to section 1.3 where in terms of the variable 


= 2sin u/2 (6) 
we have by 1(35) 
n 
(7) 
2 sin u/2 n=0 


Also by 1(39) 


v 
u/2 = arcsin v/2 = — —- 
/ / 2 23 5 


On substituting (8) into the exponential series we find the expansion 


= 2) (9) 


n=0 


with positive coefficients, the first three of which are found to be 
2 


t t 
éo(t) = |, é2(t) = e,(t) = 32 48 (10) 


On multiplying the series (7) and (9) we obtain the expansion (4). From the values 
1(40) and (10) we readily find 


| 

— 

} 

| 
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(k) (k) 


k t (k) 5k? + tk t 
d d t =—-+—, d t) = 


—+—+—. (11 
760 96 32 


Our arguments of section 2.2 may now be repeated leading to the following theo- 
rem: 
THEOREM 4. Let k be a positive integer and m an integer such that0 <2m <k+2. Then 
—1, (k) 2m—2 


Lim(x, t) = Mi(x, t) — dy Mi(x, + (— 1)” dae Mi(x) (12) 


gives rise to a “smoothing” analytic interpolation formula 


F(x) = yalem(x — n, (13) 
which is exact for the degree min(2m—1, k—1). Moreover (13) always preserves the degree 
k—1. 

The adjective “smoothing” was purposely written in quotation marks in order to 
indicate that there is no general proof as yet that (13) always reduces, for integral 
values of the variable x, to a smoothing formula in the sense of our Definition 1 of 
Part A, section 1.1. For indeed, (2) and (3), which imply such a proof, were only con- 
jectured. In the Appendix we give 8-place tables of the three basic functions 


Li(x) = L42(x, 1/8), 
Lo(x) = L42(x, 1/2), (14) 
L3(x) = Le,s(x, 1/2), 
as well as 7-place tables of their first and second derivatives. For these three sets of 
values of the parameters k, ¢, and m, the interpolation formula (13) is indeed a smooth- 


ing formula. This point is verified by an inspection of the corresponding characteristic 
functions 


o:(u) = L,(0) + 2L,(1) cos u + 2L,(2) cos 2u+---, (¢ = 1, 2, 3,). (15) 


From the values of L;(m), as given by our tables, we computed the following table 
for these characteristic functions: 


u oi(u) $2(u) ¢3(u) 
0° 1.00000 1.00000 1.00000 
30° .99734 .99519 .99952 
60° .96332 93655 .97760 
90° 85492 .76500 84693 
120° .67727 .51297 .56702 
150° 50474 . 29296 27879 
180° 43283 .20728 16123 


Since 0<¢,(u) <1 for 0<uS180°, all three formulae (13) are smoothing formulae 
according to Part A, section 1.1. Also ¢2(u) <¢:(u) implies that Le gives a stronger 
smoothing formula as compared to J. 
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Our set of tables is intended mainly for the purpose of illustrating the method. 
A more complete set of tables would be needed in order to furnish smoothing of a de- 
sired strength, as required by the needs of the numerical data at hand. 


APPENDIX 


Description of the tables and their use for the analytic approximation of equidis- 
tant data. In the Tables I, II, and III, we have tabulated the following three functions 


L(x) = M 1/8) 1/8) (1) 
L2(x) = M4(x, 1/2) — 1/2) (2) 

a(x) = Me(x, 1/2) — 1/2) + 1/2) (3) 


and their first two derivatives. The function M;(x, ¢) occurring in these definitions 
may be defined by the integral 
sin u/2\* 
t) =— e~#(u/2)" | _______ } cos uxdu. 


2a u 


A given sequence of equidistant ordinates 
{yn (4) 


is approximated by either one of the three analytic functions 


The choice among these approximations depends on the amount of smoothing de- 
sired. The formula (5), for i=1 and i=2, is exact for (i.e., reproduces) cubic polyno- 
mials. For i=3 the formula (5) is exact for quintic polynomials. For the same data (4), 
the sequence { F,(n)} is always smoother than the sequence { F:\(n)}. Generally, the 
sequence { F3(n) } should be smoother than the sequence { F,(n) }. 

The first and second derivatives of the approximation (5) may be computed by the 


similar formulae 
F{(x) = (x — »), (6) 
Fi (x) = (x — »). (7) 


The arrangement of our tables is such as to facilitate the computation of F,(x) 


by (5), as explained in the Appendix to Part A. 

An example of smoothing with subtabulation to tenths. We propose to compute a 
table of the approximation F.(x), in the range 31 <x 334, for the same ordinates 
{y,} as were used in our example of Part A (Appendix). The ordinates which we now 


require are given by the following table: 


4 
= 
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n Yn A A? A3 A‘ as 
26 32840 
27 34790 1950 
28 37260 2470 520 
29 40440 3180 710 190 
30 44750 4310 1130 420 230 
31 51120 6370 2060 930 510 280 
32 59390 8270 1900 — 160 — 1090 — 1600 
33 67550 8160 — 110 —2010 —1850 — 760 
34 73820 6270 —1890 —1780 230 2080 
35 77830 4010 —2260 — 370 1410 1180 
36 80240 2410 — 1600 660 1030 — 380 
37 || 81660 1420 — 990 610 — 50 — 1080 
38 82330 670 — 758 240 — 370 — 320 
39 || 82680 350 — 320 430 190 560 


From these values and the Table II of L(x) and L#’ (x) we obtain the following tables 
of the approximation F(x) and its second derivative Fy’ (x) shown with their differ- 
ences. 


x F(x) A | A? | A3 At Fi (x) A A? 
31.0 | 51232.76 1901.77 
31.1 | 51989.86 | 75710 1767.20 | —13457 
31.2 | 52764.62 | 77476 1766 1611.60 | —15560 | —2103 
31.3 | 53555.48 | 79086 1610 | —156 1436.84 | —17476 | —1916 187 
31.4 | 54360.69 | 80521 1435 | —175 | —19 1245.42 | —19142 | —1666 |} 250 
31.5 | 55178.34 | 81765 1244 | —191 | -16 1040.29 | —20513 | —1371 295 
31.6 | 56006.39 | 82805 1040 | —204 | —13 824.64 | —21565 | —1052 319 
31.7 | 56842.68 | 83629 824 | —216 | —12 601.66 | —22298 | — 733 319 
31.8 | 57684.98 | 84230 601 | —223| — 7 374.38 | —22728 | — 430| 303 
31.9 | 58531.03 | 84605 375 | —226| — 3 145.59 | —22879 | — 151 279 
32.0 | 59378.53 | 84750 145 | —230| — 4 || — 82.29 | —22788 91 242 
32.1 | 60225.21 | 84668 | — 82 | —227 3 || — 307.15 | —22486 302 211 
32.2 | 61068.82 | 84361 | — 307 | —225 2 || — 527.12 | —21997 489 187 
32.3 | 61907.17 | 83835 | — 526 | ~—219 6 || — 740.51 | —21339 658 169 
32.4 | 62738.12 | 83095 | — 740 | ~—214 5 || — 945.74 | —20523 816 158 
32.5 | 63559.62 | 82150 | — 945 | —205 9 || —1141.23 | —19549 974 158 
32.6 | 64369.71 | 81009 | —1141 | —196 9 || —1325.45 | —18422 1127 153 
32.7 | 65166.57 | 79686 | —1323 | —182 14 || —1496.90 | —17145 1277 150 
32.8 | 65948.46 | 78189 | —1497 | —174 8 || —1654.18 | —15728 1417 140 
32.9 | 66713.83 | 76537 | —1652 | —155 19 || —1796.14 | —14196 1532 115 
33.0 | 67461.25 | 74742 | —1795 | —143 12 || —1921.91 | —12577 1619 87 
33.1 | 68189.46 | 72821 | —1921 | —126 17 || —2030.98 | —10907 1670 51 
33.2 | 68897.38 | 70792 | —2029 | —108 18 || —2123.19 | — 9221 1686 16 
33.3 | 69584.08 | 68670 | —2122 | — 93 15 || —2198.74 | — 7555 1666 | — 20 
33.4 | 70248.81 | 66473 | —2197 | — 75 18 || —2258.04 | — 5930 1625 | — 41 
33.5 | 70890.96 | 64215 | —2258 | — 61 14 || —2301.69 | — 4365 1565 | — 60 
33.6 | 71510.12 | 61916 | —22909 | ~ 41 20 || —2330.31 | — 2862 1503 | — 62 
33.7 | 72105.98 | 59586 | —2330 | — 31 10 || —2344.57 | — 1426 1436 | — 67 
33.8 | 72678.41 | 57243 | —2343 | — 13 18 || —2345.07| — 50 1376 | — 60 
33.9 | 73227.40 | 54899 | —2344| ~ 1 12 || —2332.47 1260 1310 | — 66 


73753.07 | 52567 | —2332 12 13 || —2307.47 2500 1240 | — 70 
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A comparison of the approximation F;(x) with the strictly interpolating function 
F(x), obtained in the Appendix of Part A, is of interest. The function F(x) was ob- 
tained by the formula 

F(x) = >> yLi(x—v,), (k= 4,¢= 1/2), 


y= — 0 


where, in view of III(1) and III(2), we may define Z;,(x, t) by the expansion 


Lilx, t) = Milx, t) — Mila, t) +04 (8 Milz, — (8) 


Our present approximation F2(x) was computed by the formula (5), for i=2, where 
L2(x), by (2), I[I(5) and III(11), happens to be identical with the sum 


= Mi(x,t) — (8 Mila, 
of the first two terms of the series (8). A comparison of the tables of F(x) and F2(x) 


shows that their difference in the range 31 Sx <34 nowhere exceeds 0.23% of the 


value of F(x). 
TABLE I. L;(x) =L4,2(x, 1/8), Li (x), Li (x) 


| x+.0 x+.3 x+.4 

3 — .00041123 — .00018550 — .00007621 — .00002830 — .00000943 

2 — .03462580 — .02756299 — .02099452 — .01533104 — .01073240 

1 14220425 .08277004 03516336 — .00066601 — .02556532 

0 78566556 .77475976 . 74281594 69204028 .62577328 
14220425 .21251290 .29166041 .37661030 .46353840 
of — .03462580 — .04145730 — 04699226 — .04985257 — .04841746 
nf — .00041123 — .00083691 — .00157669 — .00277249 — .00458632 
— .00000003 — .00000017 — .00000074 

x | x+.5 | x+.6 | x+.7 x+.8 x+.9 

3 — .00000280 — .00000074 — .00000017 — .00000003 

2 — .00718751 — .00458632 — .00277249 — .00157669 — .00083691 

1 | —.04092087 — .04841746 — .04985257 — .04699226 — .04145730 

0 | .54811118 .46353840 .37661030 .29166041 .21251290 
of 54811118 62577328 69204028 .74281594 .77475976 
| — 04092087 — 02556532 — .00066601 .03516336 .08277004 
— .00718751 — .01073240 — 01533104 — .02099452 — .02756299 
—4 — .00000280 — 00000943 — 00002830 — 00007621 — .00018550 


; 
| | 
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1946] 
Li (x) 
x x+.0 x+.1 x+.2 x+.3 x+.4 
3 .0030922 .0015611 .0007155 .0002953 .0001089 
2 .0709642 .0690870 .0616092 .0514079 -0405956 
1 — .6512052 — .5359047 — .4163489 — .3016954 . 1986416 
0 .0000000 — .2168063 — .4183406 — .5915396 . 7269119 
—1 .6512052 . 7515608 . 8262912 8662852 . 8650057 
—2 — .0709642 — .0638857 — .0445045 — .0099839 .0416488 
—3 — .0030922 — .0056121 — .0094213 — .0147669 .0217952 
—4 — .0000001 — .0000006 — .0000026 .0000103 
Li (x) 
x x+.5 x+.6 x+.7 x+.8 x+.9 
3 .0000357 -0000103 .0000026 .0000006 .0000001 
2 .0304945 .0217952 .0147669 .0094213 .0056121 
1 —.1113051 — .0416488 .0099839 .0445045 .0638857 
0 — .8188067 — .8650057 — .8662852 — .8262912 . 7515608 
—1 . 8188067 . 7269119 -5915396 .4183406 . 2168063 
—2 .1113051 . 1986416 . 3016954 .4163489 .5359047 
—3 — .0304945 — .0405956 — .0514079 — .0616092 .0690870 
—4 — .0000357 — .0001089 — .0002953 — .0007155 .0015611 
Ly" (x) 
x x+.0 x+.1 x+.2 x+.3 x+.4 
4 — .0000004 — .0000001 
3 — .0197354 — .0114013 — .0059477 — .0027759 .0011500 
2 .0202421 — .0521077 — .0926044 — .1078983 . 1061742 
1 1.0966561 1.1912649 1.1844282 1.0974670 .9568292 
0 —2.1943249 —2.1159568 —1.8927238 —1.5554043 . 1427100 
—1 1.0966561 . 8923803 . 5868711 . 2020705 . 2337191 
—2 .0202421 . 1269947 . 2654156 .4283134 -6058876 
—3 — .0197354 — .0311722 — .0454298 — .0617347 .0788289 
—4 — .0000004 — .0000019 — .0000092 — .0000377 .0001346 
x x+.5 x+.6 x+.7 x+.8 x+.9 
3 — .0004201 — .0001346 — .0000377 — .0000092 .0000019 
2 — .0947351 — .0788289 — .0617347 — .0454298 .0311722 
1 . 7867260 .6058876 .4283134 . 2654156 . 1269947 
0 — .6915708 — .2337191 . 2020705 . 5868711 . 8923803 
— .6915708 —1.1427100 —1.5554043 —1.8927238 . 1159568 
—2 . 7867260 - 9568292 1.0974670 1.1844282 . 1912649 
—3 — .0947351 — .1061742 — .1078983 — .0926044 .0521077 
-—4 — .0004201 — .0011500 — .0027759 — .0059477 .0114013 
—5 .0000001 
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TABLE II. Lo(x) =L42(x, 1/2), Lz (x), La’ (x) 
L2(x) 
x x+.0 x+.1 x+.2 x+.3 x+.4 

5 |  —.00000001 

4 | —.00003297 — .00001721 — .00000871 — 00000428 .00000203 

3 — 00453670 — 00313034 — 00210474 — 00137859 .00087930 

2 — .04033977 — .03775949 — .03379714 — .02913417 02429574 

1 20271717 . 14753253 09914645 .05829935 .02523450 

0 68438455 67711287 65567464 62117134 57534514 
—1 .20271717 . 26343912 .32793367 . 39399436 45907812 
—2 — .04033977 — .04070831 — .03791256 — .03092183 01869149 
—3 — .00453670 — .00640787 — .00882100 — .01183236 .01545906 
—4 — 00003297 — 00006127 — .00011052 — 00019364 .00032972 
-5 | —.00000001 — .00000003 — .00000007 — .00000018 00000042 

x | x+.5 x+.6 x+.8 x+.9 

4 | —.00000094 — .00000042 — .00000018 — .00000007 .00000003 

| —.00054590 — 00032972 — 00019364 — .00011052 .00006127 

2 — 01965692 — .01545906 — .01183236 — .00882100 .00640787 

1 — .00024449 — 01869149 — 03092183 — .03791256 .04070831 

0 . 52044824 45907812 . 39399436 . 32793367 26343912 
. 52044824 57534514 62117134 65567464 67711287 
— 00024449 .02523450 .05829935 09914645 .14753253 
—3 — .01965692 — .02429574 — .02913417 — .03379714 03775949 
—4 — 00054590 — .00087930 — .00137859 — .00210474 .00313034 
—5 — .00000094 — .00000203 — .00000428 — .00000871 00001721 

Li (x) 

x | x+.0 | x+.2 x+.3 x+.4 

s | .0000001 | 

4 | .0002093 0001145 .0000607 .0000311 0000154 

0162494 0120195 | .0086291 0060153 .0040721 

.0162409 .0339777 0441321 0482532 .0478931 

1 | —.5820676 —.5195203 | —.4469715 — .3695748 . 2920635 

0 | .0000000 —.1448007 | —.2821139 — .4050264 .5077160 
-1 | . 5820676 6294198 | 6567760 6601752 .6369100 
—2 | —.0162409 0104651 | 0471784 .0943909 . 1518602 
—3 | —.0162494 —.0213049 | —.0270545 — 0332049 .0392596 
—4 | —.0002093 — .0003705 .0006358 — .0010581 .0017083 
—5 | —.0000001 —.0000002 | —.0000006 — .0000015 .0000034 


= 
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Li (x) 

x+.5 x+.6 x+.7 x+.8 x+.9 
.0000074 .0000034 -0000015 .0000006 .0000002 
.0026769 .0017083 .0010581 .0006358 .0003705 
.0444848 .0392596 .0332049 .0270545 .0213049 
. 1876748 — .1518602 — .0943909 — .0471784 — .0104651 
. 5858619 — .6369100 — .6601752 — .6567760 — .6294198 
5858619 .5077160 -4050264 . 2821139 . 1448007 
. 1876748 . 2920635 . 3695748 -4469715 .5195203 
.0444848 — .0478931 — .0482532 — .0441321 — .0339777 
.0026769 — .0040721 — .0060153 — .0086291 — .0120195 
.0000074 — .0000154 — .0000311 — .0000607 — .0001145 

Lz" (x) 
x+.0 x+.1 x+.2 x+.3 x+.4 
.0000010 — .0000004 — .0000001 

.0012284 — .0007087 — .0003952 — .0002129 — .0001108 
.0465343 — .0380479 — .0298698 — .0225871 — .0164844 
. 2201282 1369510 .0687454 .0162834 — .0210828 
.5579563 .6842412 . 7580806 - 7819152 . 7615818 
1.4606815 —1.4227560 —1.3118991 —1.1365778 — .9099683 
.5579763 3809905 . 1594624 — .0960473 -3711575 
. 2201282 -3157558 -4193623 .5245097 . 6230076 
.0465343 — .0543645 — .0601451 — .0620423 — .0578316 
.0012284 — .0020585 — .0033357 — .0052283 — .0079270 
.0000010 — .0000025 — .0000057 — .0000126 — .0000270 

x+.5 x+.6 x+.7 x+.8 x+.9 
.0000557 — .0000270 — .0000126 — .0000057 — .0000025 
.0116249 — .0079270 — .0052283 — .0033357 — .0020585 
.0450247 — .0578316 — .0620423 — .0601451 — .0543645 
7053804 .6230076 .5245097 -4193623 .3157558 
.6486751 — .3711575 — .0960473 . 1594624 . 3809905 
.6486751 — .9099683 —1.1365778 —1.3118991 —1.4227560 
. 7053804 . 7615818 . 7819152 . 7580806  ,6842412 
.0450247 — .0210828 .0162834 .0687454 . 1369510 
0116249 — .0164844 — .0225871 — .0298698 — .0380479 
.0000557 — .0001108 — .0002129 — .0003952 — .0007087 
— .0000001 — .0000004 
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TABLE III. Ls(x) =Le,3(x, 1/2), Ls (x), La’ (x) 
x | x+.0 x+.1 2+.3 x+.4 
6 | .00000024 .00000012 .00000005 .00000002 .00000001 
5 .00008286 .00005169 00003150 .00001874 .00001088 
4 | .00268003 .00209760 .00160285 .00119677 .00087361 
3 .00146844 .00384345 .00530620 .00602933 .00618622 
2 .06657958 .06134625 — .05394260 — .04532872 — .03631488 
1 | . 20814167 . 14737265 09315196 -04655844 .00818720 
0 . 70841267 . 70114512 -67967926 .64500315 .59869152 
—1 . 20814167 . 27390708 . 34268907 .41215975 .47975381 
—2 .06657958 .06856746 — .06617125 — .05825575 — .04376882 
—3 .00146844 .00197317 — .00659556 — .01244745 — .01948472 
—4 .00268003 .00333890 .00404886 .00476644 .00542627 
—5 .00008286 .00012978 .00019870 00029745 .00043550 
—6 .00000024 .00000049 . 00000096 .00000183 .00000340 
| x+.5 x+.6 x+.7 x+.8 x+.9 
5 .00000616 .00000340 .00000183 .00000096 .00000049 
4 .00062367 .00043550 00029745 .00019870 .00012978 
3 .00593823 .00542627 .00476644 .00404886 .00333890 
2 .02754330 .01948472 — .01244745 — .00659556 — .00197317 
1 .02182865 .04376882 — .05825575 — .06617125 — .06856746 
0 54280388 -47975381 -41215975 34268907 .27390708 
=f 54280388 59869152 64500315 .67967926 .70114512 
02182865 .00818720 04655844 09315196 14737265 
02754330 03631488 — .04532872 — .05394260 — 06134625 
.00593823 00618622 .00602933 00530620 00384345 
~§ 00062367 00087361 .00119677 00160285 .00209760 
~6 00000616 .00001088 00001874 00003150 00005169 
.00000001 .00000002 .00000005 .00000012 
Lj (x) 
x+.0 x+.1 x+.2 x+.3 x+.4 
6 .0000017 .0000009 — .0000004 — .0000002 — .0000001 
5 .0003813 .0002500 — .0001598 — .0000996 — .0000605 
4 .0062356 .0053949 — .0044992 — .0036324 — .0028472 
3 .0288514 .0189098 .0106360 .0041176 — .0007095 
.0379708 .0648953 .0815717 -0893524 .0898535 
1 .6356363 -5771501 — .5054638 — .4254184 — .3417929 
0 .0000000 . 1447858 — .2828738 — .4080067 — .5147714 
—1 6356363 -6763588 -6953825 6897327 -6576776 
—2 .0379708 .0001221 .0497076 . 1103706 1808527 
—3 .0288514 .0401754 — .0523583 — .0646105 — .0758656 
—4 .0062356 .0069002 .0072273 0070158 .0060307 
—5 0003813 .0005675 .0008239 .0011670 .0016121 
—6 .0000017 .0000033 . 0000063 -0000116 .0000206 
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Li (x) 
x x+.5 x+.6 x+.7 x+.8 x+.9 
5 — .0000358 — .0000206 — .0000116 — .0000063 — .0000033 
4 — .0021709 — .0016121 — .0011670 — .0008239 — .0005675 
3 — .0040148 — .0060307 — .0070158 — .0072273 — .0069002 
2 -0847953 .0758656 -0646105 .0523583 .0401754 
1 — .2590022 — .1808527 — .1103706 — .0497076 — .0001221 
0 — .§989336 — .6576776 — .6897327 — .6953825 — .6763588 
.5989336 .5147714 -4080067 . 2828738 . 1447858 
—2 2590022 .3417929 -4254184 .5054638 -5771501 
—3 — .0847953 — .0898535 — .0893524 — .0815717 — .0648953 
—4 .0040148 .0007095 — .0041176 — .0106360 — .0189098 
.0021709 .0028472 -0036324 -0044992 .0053949 
—6 .0000358 .0000605 .0000996 .0001598 .0002500 
.0000001 .0000002 .0000004 .0000009 
(x) 
x x+.0 x+.1 x+.2 x+.3 x+.4 
6 .0000117 .0000061 -0000031 .0000015 .0000007 
5 .0015631 -0010868 0007350 . .0004835 .0003094 
4 .0077612 .0088518 .0089242 .0083240 .0073357 
3 — .1069900 — .0913730 — .0739648 — .0565175 — .0403140 
2 . 3239975 -2160536 . 1197566 .0385786 — .0256675 
1 .5032414 -6588169 7667843 .8261081 8390021 
0 —1.4591699 —1.4253081 —1.3259553 —1.1676245 — .9605716 
-1 .5032414 . 3045939 .0708748 — .1867923 — .4548011 
—2 . 3239975 -4384255 -5526106 .6585611 . 7474910 
—3 — .1069900 — .1186224 — .1237228 — .1195410 — .1033391 
—4 .0077612 .0052612 .0009447 — .0055676 — .0145672 
—5 .0015631 .0021860 :0029705 .0039174 .0050049 
—6 .0000117 .0000217 _.0000392 .0000685 .0001165 
—7 .0000001 .0000001 
x x+.5 x+.6 x+.7 x+.8 x+.9 
6 .0000003 .0000001 .0000001 
5 .0001926 .0001165 .0000685 .0000392 .0000217 
4 .0061776 .0050049 .0039174 .0029705 .0021860 
3 — .0261859 — .0145672 — .0055676 .0009447 .0052612 
2 — .0726654 — .1033391 — .1195410 — .1237228 — .1186224 
1 .8104398 7474910 .6585611 .5526106 .4384255 
0 — .7179590 — .4548011 — .1867923 — .0708748 3045939 
-1 — .7179590 — .9605716 —1.1676245 —1.3259553 —1.4253081 
—2 . 8104398 . 8390021 .8261081 . 7667843 -6588169 
—3 — .0726654 — .0256675 .0385786 .1197566 . 2160536 
—4 — .0261859 — .0403140 — .0565175 — .0739648 — .0913730 
—5 .0061776 .0073357 .0083240 .0089242 .0088518 
—6 .0001926 .0003094 .0004835 .0007350 .0010868 
.0000003 .0000007 .0000015 .0000031 .0000061 
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AN ITERATION METHOD FOR CALCULATION 
WITH LAURENT SERIES* 


BY 


H. A. RADEMACHER anp I. J. SCHOENBERG 
University of Pennsylvania and Ballistic Research Laboratories, Aberdeen Proving Ground 


Introduction. The power series is a basic concept of Analysis which is of funda- 
mental importance from the theoretical as well as from the computational point of 
view. The theoretical importance of power series springs from the fact that it repre- 
sents any analytic function in the neighborhood of a regular point. The reason for its 
practical importance is the ease with which implicitly defined functions, by finite 
relations or differential equations, may be expanded in power series by the so-called 
method of undetermined coefficients, known and used since the dawn of mathematical 
analysis. 

Laurent series play a definitely minor role as compared to power series. One reason 
is the more complicated nature of the connection between the sum of the series and 
its coefficients. Another reason, dependent on the first, is the difficulty of calculations 
with Laurent series. 

The purpose of this paper is to describe a method whereby rational or algebraic 
operations with Laurent series may be performed with high accuracy at the expense 
of a reasonable amount of labor. A general approximation method to empirical data, 
developed by one of us,’ required the very accurate reciprocation of certain Laurent 
series. This problem of reciprocation of Laurent series was the starting point of our 
investigation. Our method for solving this particular problem turned out to be identi- 
cal with a method of reciprocation of finite matrices already investigated by 
H. Hotelling.? We finally point out that our method of computation with Laurent 
series extends to computations with trigonometric series provided these series con- 
verge absolutely. 

1. Newton’s algorithm and statement of the problem. Let 


f(x) = dox™ of +--+ + an = 0, (do 0), (1) 


be an algebraic equation with numerical real or complex coefficients. If x is a simple 
root of this equation, then very close approximations to x may be readily computed 
by Newton's iterative algorithm represented by the recurrence relation 


S(%n) 


The reason for the fast convergence of x, towards x is as follows: Expanding the right- 


(2) 


Xn+1 = Xn 


* Received Nov. 6, 1945. 

1], J. Schoenberg, Contributions to the problem of approximation of equidistant data by analytic func- 
tions, Part A, Quart. Appl. Math. 

2H. Hotelling, Some new methods in matrix calculation, Ann. Math. Statist. 14, 1-34 (1943), espe- 
cially p. 14, and Further points in matrix calculation and simultaneous equations, Ann. Math. Statist. 14, 


440-441 (1943). 
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hand side in a power series in x, —x, in the neighborhood of the simple root x, we find 
that (2) may be written as 


Xng1 — X = Co(Xn — + — x)P+---, (3) 


with coefficients c, depending only on the root x. Because there is no linear term 
in x,—x on the right-hand side we find that from a certain point on the error X,41—* 
is of the order of magnitude of the square of the previous error x,—x. From this stage 
on, each step will approximately double the number of correct decimal places of the 
previous approximation x,. This type of rapid convergence is sometimes referred to 
as “quadratic convergence.” 

Let us notice that the iterative process (2) requires a division, by f’(x,), at each 
step of the process. This is a serious handicap in computing with machines which do 
not perform the operation of division, as for example the standard punch-card ma- 
chines. This division is likewise a handicap if we wish to extend the process to the 
realm of matrices where division is a difficult numerical operation. 

We propose to modify Newton’s algorithm (2) so as to require only the operations 
of addition, subtraction and multiplication in its performance. It will then be shown 
how the modified Newton algorithm allows us to carry out numerically rational as 
well as algebraic operations on Laurent series. The most general numerical problem 
whose solution is facilitated by our method may be formulated as follows. 


PROBLEM. Let 
2) = ao(z)w™ + +--+ + an(z) = 0 (4) 


be an equation with the following properties: 
1. The coefficients a,(z) are all regular and uniform functions of 2 in the ring 


R: n<|2| (5) 
2. We have 
a(z) #0 in R. (6) 
3. The discriminant D(z) of (4) satisfies 
D(z) #0 in R (7) 


so that the equation (4) has no critical point in R. Let now w=w(z) be a branch of a solu- 
tion of (4) which is necessarily regular in R but need not be uniform in R. Given the nu- 
merical values of the Laurent expansions of the coefficients a,(z), the problem is to find the 
values of the coefficients of the Laurent expansions of w(z). 


Remarks. 1. The difficulty of this problem is due to its being concerned with 
Laurent series rather than ordinary power series. Indeed, if everything else is un- 
changed, we replace, in its formulation, the ring (5) by the circle | z| <re, then all 
Laurent series mentioned become power series, in which case the power series expan- 
sion of the branch w=w(z) may be obtained by the method of undetermined coefh- 
cients (see the first paragraph of our Introduction). 

2. We did not require the branch w=w(z) to be uniform in R. However, we do 
not restrict our problem by assuming w(z) to be uniform in R. Indeed, if w(z) returns 
to its values after k turns in R, k>1, we change variable by setting 
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Our equation (4) thereby becomes 
f(w, =0 


and the branch w(z) becomes uniform in the corresponding ring in the ¢-plane. If we 
can determine its uniform Laurent series 


w(z) = wnt” 
we also have its algebraic expansion 


ll 
€ 
= 
v2 

2 
= 


w(z) 


3. Even the case m=1 is far from trivial. Thus 
a(z)w—-1=0 


amounts to the important problem of the reciprocation of a given Laurent series. 
2. The modification of Newton’s algorithm. We return in this section to the case 
of the ordinary algebraic equation (1). We now impose the restriction that 


f(x) has only simple zeros. (8) 


This condition implies that the polynomials f(x), f’(x) have no common divisors and 
that we can therefore determine uniquely, by rational operations alone, two poly- 
nomials $(x) and W(x) satisfying the identity 


S(x)o(x) + f'(x)¥(x) = 1, (9) 


and such that the degrees of ¢ and w do not exceed m—2 and m—1, respectively. 
The coefficients of ¢(x), ¥(x) are rational functions of the coefficients a,. For later 
reference it is important to remark that the coefficients of ¥(x) may be written asa 
quotient of polynomials in a, divided by the discriminant D of the polynomial f(x). 
Indeed, the coefficients of ¢ and y are determined, in view of (9), by a system of linear 
equations whose determinant is precisely the discriminant D of f(x). This procedure 
leads to explicit expressions of y and D in determinant form. Thus for m =3 we obtain 


| 1 0 0 | 3a 0 a 
‘ | 3a 0 O a O | 2a, 3a O a 4% 
¥(x) = > | 2a, 3a O a, ai, D = 2a; ade (10) 
| O 2a, as ae! 


This expression, which generalizes to any value of m, indeed shows that the coeffi- 
cients of ¥(x) have the common denominator D if regarded as rational functions of 


the a,’s. 
Now we modify Newton’s algorithm (2) to its new form’ 


3 We learn from a note by J. S. Frame, Remarks on a variation of Newton’s method, Amer. Math. 
Monthly, 52, 212-214 (1945), that precisely the same modification of Newton’s algorithm has already 
been used since 1942 by H. Schwerdtfeger, of the University of Adelaide, South Australia, for the numeri- 
cal solution of ordinary algebraic and transcendental equations. 


| 
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= Xn — S(%n)W(%n). (11) 
Setting 
F(x) = x — f(x)¥(x) (12) 
we may write (11) as 
= F(x). (11’) 


On comparing (2) and its modification (11) we see that the division required by (2), 
at each step of the process, is not present in (11). Now we want to show that the 
algorithm (11%) also enjoys the property of (2) of producing fast convergence towards 
the zeros of f(x). Indeed, let x be a root of (1), 


f(x) = 0, (13) 
and let us expand F(x,) about the point x, =x. Writing for convenience f(x) =f, 
Y (x) =~, we have by Taylor’s formula 

f(%n) =f + — + — 
=P (an — x) + — 


hence 
By (9) and (13) we have f=0, f’f) =1 and therefore 
— = — 3(2f'V + — x)? 
This shows that we may write our relation (11) in the form 
— = bo(x)(%_ — x)? + ba(x) — +--+ + — (14) 


where the b,(x) are polynomials in x with coefficients which are polynomials in a, divided 
by the common denominator D. 

Again the missing linear term in x,—<x, on the right-hand side of (14), shows that 
if xo is sufficiently close to x, then the algorithm (11) will insure that x,—x with quad- 
ratic convergence. 

An important special case of (1) is the equation 


ax™—1=0, (a ~ 0). (15) 
The identity 
x 
(— 1)\(ax™ — 1) + — = 1 
m 


shows that in this case 


¥(x) = — x. 
m 


The relation (11) now becomes 
1 
= Xn + — — axe). (16) 
m 


In particular if m=1, (15) reduces to 
ax—1=0 (17) 


% 
aa 
= 
7 
a 
| 
4 
| 
1 
an 
j 
| 
| 
4, 
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when (16) becomes 
= Xn + x,(1 — (18) 

3. The reciprocation of matrices. The advantage of the modified, or division-free, 
Newton algorithm (11) appears in connection with matrix calculations. In recent 
years H. Hotelling has recommended the following procedure of finding the reciprocal 
X =A7 of a given numerical non-singular matrix 


A=llaj| =1,--+,m). (19) 


Obtain in some way, e.g. by the so-called Gauss, or Doolittle, process a good approxi- 
mation X» to A~!. Then improve this approximation by the recurrence relation 


= Xn + — AX,). (20) 


In the case of m =1 this relation is identical with (18). 
In studying the convergence of X, towards X = A— Hotelling metrizes the space 


of real m Xm matrices by means of the absolute value or norm 


N(A) = 4/ Das (21) 


which enjoys the following properties 
N(A + B) S N(A) + N(B), = N(AB) S N(A)N(B). (22) 
By means of these inequalities Hotelling derived an estimate of N(X,—X) which 
was improved by A. T. Lonseth as follows: 
INEQUALITY OF HOTELLING AND LONSETH. Let Xo be an approximation to X =A 
such that 


NU — AX) =k <1. (23) 


Starting with Xo we obtain the sequence X,, by (20). Then 
N(X, — X) S N(Xo)-k™-(1 — (24) 


This interesting result shows in particular that the inequality (23) is sufficient to 
insure the convergence of the process. 

Our generalizations (16) and (11) of the recurrence relation (18) suggest similar 
iterative procedures for the solution of non-linear algebraic matrix equations. We 
prefer, however, to pass on to a discussion of calculations with Laurent series. 

4. Calculations with Laurent series. Let 


a(z) = an2", (n< | 2| < 12), (25) 


be a Laurent series converging in the ring (5). There is no inherent restriction of the 
generality of the Problem formulated in our Introduction if we assume that the-ring R 
contains the unit-circle |z| =1, i.e. 


4 See Hotelling’s second note already mentioned. 


: 
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<1 (26) 


An advantage of this normalization is that it implies that a,—0 exponentially as 
n—-+ © or n—+— ©, insuring that the sequence {a,} is “finite” to a fixed number of 


decimal places. 
The relation (25) sets up a one-one correspondence 


a(z) ~ {an} 


between functions a(z) uniform and regular in R and sequences {a,} for which the 
series (25) converges in R. To the function a(z) =1 corresponds the unit-sequence 


I: ao =1, a =0 if n¥0. 


This correspondence may be interpreted as an isomorphism concerning the operations 
of addition, subtraction, multiplication and multiplication by a scalar. Indeed, if 


b(z) = (7, <2 < 12), (27) 


is a second series then we find, on multiplying (25) and (27) that to the product 


c(z) = a(z)b(z) (28) 
corresponds the series 
c(z) = (29) 
where 
Ta On—vBy. (30) 


Thus to the operation (28) of multiplication of the functions a, b, corresponds the 
operation of convolution (30) of the two sequences {a,}, {8,}, an operation which 
we write as 


y = af. (31) 


We mention incidentally a third interpretation of Laurent series isomorphic to the 
two already discussed. Indeed, consider the (4-way) infinite matrix 


(32) 


in which a;_; is the element in the ith row and jth column, both 7 and j assuming 
all integral values. Such matrices may be designated “striped” for the reason that all 
elements lying on a line, sloping down at a 45° angle, are identical. To every sequence 
{a,} corresponds one such matrix, and conversely. The isomorphism between such 
matrices and sequences becomes evident if we remark that the multiplication of two 
striped matrices ||a;~;||, is another striped matrix where the sequence 
{yn} is given by (30). This remark throws some light on the connection between 
Laurent series and the case of finite matrices of §3. 


4 
7 
‘ 
a 
= 
— 
Ad 
7 
4 
“4 
7 
| 
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We finally define the norm of the function a(z), or of the sequence {a,} , as the 
non-negative number 


= N(a) = (33) 


This norm also enjoys the two properties (22) or 
N(a + b) S N(a) + N(d), N(ab) S N(a)N(b). (34) 


Their verification is immediate in this case. 

We are now able to attack the general Problem of section 1. However, it is essen- 
tial to discuss first the important case m =1 of reciprocation. 

4.1. The reciprocation of Laurent series. With the norm of a Laurent series as 
defined by (33), the result of Hotelling and Lonseth (section 3) applies to Laurent 
series without any change. Assuming that the sum a(z) of the given Laurent series 
(25) does not vanish in R, we are to find the expansion 


1 
w(z) = = (35) 
a(z) 
Let 
Wo(z) = as" (36) 


be an approximation to (35) such that 


N(i — aw) = NU — aw) = k <1. (37) 


The very important problem of how such approximations may be obtained will be 
discussed later (section 4.3). This starting sequence is now to be improved by the 


relation 


The rapid convergence is assured by the Hotelling-Lonseth inequality 
N(w™ — w) S N(w)-k*-(1 — Rk). (39) 


Pending a discussion of procedures for obtaining the first approximation, we may 
therefore regard the numerical problem of reciprocation as solved. This implies that 
we may perform all four rational operations on Laurent series and that we may thus 
find the Laurent expansion of any rational function of Laurent series. 

4.2. The general algebraic case. We turn now to the general case of the equation 
(4) with the two additianal, and as we have seen, unessential restrictions that our 
ring R contains the unit-circle | z| =1 and that the solution w=w(z), of (4), be uni- 
form in R. The problem is to find the numerical values of the coefficients of the 


expansion 


w(z) = (40) 


We return to our discussion (section 2) of the division-free Newton algorithm (11), 
especially in its expanded form (14). This discussion remains valid if applied to (4) 
rather than (1). The algorithm is in this case 
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Watt = Wn — S( Wa, z)W(Wn), (41) 

the expanded form of which is 
Wnt1 — W = bo(w)(We — w)? + +++ + — (42) 


We have to remember, however, that ¥(w,) is a polynomial in w, with coefficients 
which are polynomials in a, divided by D=D(z). Since D is a polynomial in the a,, 
we may first derive its Laurent expansion by additions and multiplications from the 
given Laurent series of the coefficients a,(z) of (4). Secondly, since 


D(z) #0 in R, 


we may also find by the method of section 4.1 the Laurent expansion of 1/D(z). In 
this way we arrive at the Laurent expansions of the coefficients of ¥(w,). These pre- 
liminary Laurent series operations allow to put the relation (41) in the form 


2m—1 


Wn41 = Wa + (co(z) + €1(Z) Wn + + Com—1(Z) Wn ), (43) 
where 
are numerically known Laurent expansions. 
Now let 
wi(s)= s (45) 


be an approximation to (40), (See section 4.3.) Starting from this approximation we 
obtain the successive series 

n)v 
walt) = Dw," 2 (46) 


by means of (43). This operation of deriving w,41 from w, is of course to be performed 
on the corresponding sequences of coefficients. By (43), (44), (46), the operation takes 
the form 


Will the expansion (46) converge towards the expansion (40) of the solution? To 
answer this question we return to the form (42) of our relation. Taking the norms of 
both sides of (42) and using the properties (34) of the norm, we obtain 


— w) S N(bo(w))[N(we — w)]? +--+ + [N (wn — (48) 


This relation shows that if 
— w) = >> @, | (49) 


is sufficiently small then (48) will indeed imply 
lim N(w, — w) = 0 (50) 


no 


with quadratic convergence. 


q 
— 
a 
4 
a 
| 
| 
| a 
= 
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4.3. Derivation of an approximate Laurent expansion. The method for computa- 
tion with Laurent series described in the previous sections will now become effective 


provided we can solve the following problem. 


INITIAL APPROXIMATION PROBLEM. Let 
F(z) = >> (51) 


be regular in the ring R containing the circle | z| =1. This function F(z), whose Laurent 
coefficients c, are unknown, is defined by an algebraic equation which allows us to com- 
pute the value of F(z) for any given z of R, in particular for any root of unity. We are to 
describe a practical method whereby, given e>0, we may compute the coefficients c,* of 


a Laurent series 
F*(z) = >> (52) 


regular in R, such that 
(53) 


N(F — F*) = 


We shall now solve this problem by the method of trigonometric interpolation.5 
Let m be a positive integer and let 


Zu = (u 0, 1, 1), (54) 


be the mth roots of unity. These roots of unity satisfy the following orthogonality 


relations 
1 if v = s (mod m) 
— = 
0 if »#-s(modm). (55) 
If m is odd, m=2n-+1, we consider the Laurent polynomial 
F(z) = Cm,2” (56) 


having m arbitrary coefficients. 
If m is even, m =2n, we define our polynomial so as to contain again only m arbi- 


trary coefficients as 


n—1 


vez—(n—1) 


Whether m is even or odd we may always write 


5 Concerning the subject of trigonometric interpolation we refer to the classical memoir by Ch. J. de 
la Vallée Poussin, Sur la convergence des formules d’interpolation entre ordonées equidistantes, Bulletin de 
l’Academie royale de Belgique, 319-410 (1908), and to Dunham Jackson, The theory of approximation, 
American Mathematical Society Colloquium Publications, vol. 11, New York, 1930, chap. IV. 


J 
| 
F,,(z) = n= (58) 
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where the summation symbol },’ is to indicate that if m is even, then 


(59) 


Cm,—n = Cm,ny 


and that the terms of (58) for y= +” are to be taken with half their value. The rela- 
tion n= |m/2] is to indicate that n is the greatest integer not exceeding m/2. 

We shall now require the Laurent polynomial (58) to interpolate the function (51) 
in the points (54). This gives the m equations 


= F(z,), (u 0, 1, 1). (60) 


=—n 


On multiplying (60) by 2,/m (s fixed, —n Ss Sn) we find in all cases, in view of (55) 
after summation by yp, that 


Cae = 1 (—-nSsZn). (61) 


The construction of our approximate Laurent expansion (58), i.e. (52), has now been 
completed. The following theorem will now show that the condition (53) may also be 
realized by the present method of construction. 


THEOREM. We assume the Laurent series 


F(z) = > oz" (62) 
to converge absolutely on the unit circle |z| =1, i.e. 


< w,6 (63) 


Then our interpolating Laurent polynomial (58) satisfies the condition 


lim N(F — Fn) = 0. (64) 


mo 


Remark. Notice that the regularity of F(x) in a ring containing | z| =1 implies 
our condition (63) but not conversely. This remark is of importance concerning cal- 
culations with absolutely convergent Fourier series. (See section 5.) 

Proof. Let N be a positive integer. We shall restrict ourselves to values of m22N, 
hence »=N. We may then write 


N-1 
v=—N+1 


We shall now estimate the three sums on the right-hand side. 


6 See Dunham Jackson, loc, cit., for other conditions insuring the convergence of trigonometric 
interpolation. 


aos 
4 
— 
¥ 
} 
5 
a 
= 
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Let ¢>0 be given. In view of our assumption (63) we may choose N such that 


> (66) 


|¥|2N 


An upper bound for the second sum of (65) may now be obtained as follows. By (61) 
and (62) we have 


Cme = — >. By = — 
m v=—00 m u=0 
and finally by (55) 
Cm,s = 2. Cy (67) 
v=s(mod m) 


For all m=2N, whether m is even or odd, we now have 


Ns|s|Sn NsS|s|Sn v=s(mod m) 


Cy = € (68) 


by (66). We now return to (61). Since (62) converges uniformly on | z| =1, F(z) is 
continuous on |z| =1 and (62) is its Fourier series. We therefore have the Fourier- 
Cauchy relations 


1 
(69) 
2ri Jz|=1 


From the definition of this integral as a limit of Cauchy sums we may now write 
(with z, =1) 


1 m—-1 
¢ =lim — F(2,) 25°" — 
m— Tt p=0 


c, = lim >> 1), 


m0 poo 


and finally, by (61), 


1 m—1 
¢ = lim — >> F(z,)z,” = lim Cm». (70) 
mo M m— 2 
We now have indeed 
N-1 
<€, provided m > mo(c). (71) 
y=—N+1 


By (65), (66), (68), and (71) we now have 
N(F — Fn) < 3, provided m > mo(e), (72) 


and our theorem is established. 
4.31. The 24-ordinate scheme of numerical harmonic analysis. The interpolation 
of our given function F(z) in the 24th roots of unity will provide satisfactory approxi- 
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mation for most ordinary purposes. Let us assume for definiteness that F(z) is real for 
real z. On the unit-circle z =e” we then have 


F(e"*) = + (73) 


where the real part R(@) is an even function, while the imaginary part J(@) is odd. 
Denote by 
F, = R, + i,, (u = 0, 1,---, 23), (74) 


the computed values of our function at 15°—intervals in @, i.e. for the points (54) 
with m=24. We now interpolate the 13 ordinates R, (u=0, 1, ---, 12) by a cosine 
polynomial 


Ay + A, cos 6 +--+ + Ay, cos 116 + Aj2 cos 128, (75) 
and the 13 ordinates Jy>=0, Ih, - - , Iu, by a sine polynomial 
B,sin@ +--+ + By sin 116. (76) 


These polynomials are readily obtained by the 24-ordinate scheme as described in 
E. T. Whittaker and G. Robinson, The calculus of observations, ed. 3, 1940, section 137, 
pp. 273-278. The complex function (73) is now interpolated in the 24 points by the 
trigonometric polynomial 


Fa = Ay + A, cos 6 + +++ + Ay cos 110 + Aj2 cos 120 


+ iB; sin@ +--+ + iBy sin 110. 
Setting 
z= cos vO = + i sin v0 = 3(2” — (77) 


we obtain the Laurent sum with real coefficients 


Fo4(z) => Ao 4+ 3(A, B,)z” + 3A 


+ > 4(A, — + $A (78) 


This initial approximate Laurent expansion will be used in section 6 in our example 
of reciprocation of a Laurent series. 

5. Calculation with Fourier series. The method of calculation with Laurent series 
described in sections 4, 4.1, 4.2, 4.3 and 4.31, applies unchanged to the realm of abso- 
lutely convergent Fourier series written in the complex form 


F(z) = where z 


with the definition of the norm as 


NF) = 


The general problem of section 1 may now be reformulated, replacing the ring R by 
the unit circle | z| =1. The coefficients a,(z) of the equation (4) are now defined by 


| 
q 
q 
. 
and 
4 
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given absolutely convergent Fourier series. The conditions (6) and (7) remain un- 
changed. The important fact that a uniform, continuous solution w=w(z) of (4) 
along the unit-circle admits of an absolutely convergent Fourier expansion is now as- 
sured by a general theorem of N. Wiener and P. Lévy.’ The effectiveness of the inter- 
polation method of section 4.3 for obtaining a satisfactory initial approximate Fourier 
series is secured by our theorem of section 4.3. 

We finally mention briefly the special problem of the reciprocation of a non- 


vanishing absolutely convergent Fourier series 
= 4a0 + >> (a, cos nO + 5b, sin (79) 
n=1 


with real coefficients a,, b,. In applying our method, we have to pass to the complex 
variable z, by means of the relations (77), obtaining the series 


A(0) = a(z) = anz”, (80) 
which is now to be reciprocated. The coefficients a, being complex, it would appear 


that computations with complex numbers are unavoidable. This, however, is not the 
case since we may proceed as follows. Working with the real series 


f(z) = = Main)", (81) 
we have by (80) 
1 1 g 
= = = 82 


Starting with (81) and operating with real series only, we now form the expansions 
of f?, g? and then f?+g?. The real “Laurent” series of f?+-g? is now reciprocated by the 
method of section 4.1 and finally the series for 


HF +8), 
obtained by multiplications. This furnishes the complex values of the w, of (82). 


Returning to the variable 0, by (77) we finally obtain the ordinary real Fourier ex- 


pansion of 


1/A (6) 


6. An example of reciprocation of a Laurent series. Our numerical example will 
benefit by the following general remark concerning the modified Newton algorithm 
(11). For simplicity we limit ourselves to the case of the equation (15) or 

ax™—-1=0 (83) 
which is solved by the recurrent relation 


7 See Antoni Zygmund, Trigonometrical series, Warszawa-Lwéw, 1935, pp. 140-142. 
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1 

Xn41 = + -- — (84) 
m 


Let us assume that our first approximation x9 is of such accuracy that x2 will have all 
the accuracy we want, while x; does not quite do. More precisely we assume the 


“residual” 
(85) 


= 1 — ax,” 


so small that we may neglect r* everywhere in our calculations. We may use this fact 
in eliminating x; between the two equations 


1 

xo + — xo(1 xd"), 

(86) 


1 
m 


Indeed, by (85), (86), we have 


1 
1 = +- 
m 


m—1 
x" = +r+ r). 


and neglecting r? we find 


If we then compute x2 in this way, i.e., neglecting r* wherever it appears, we easily find 
the following approximation to x2: 


1 m+ 1 
= r). (87) 
m 


2m 


We may interpret both equations (85), (87) as a recurrence relation furnishing x? in 
terms of the first approximation x9. This process converges “cubically.” Indeed, a 
simple calculation will show that we may write (87) as 

a (m + 1)(2m + 1) 


6x? 


xf — 


(x9 — x)* + (terms of order > 3). 


We note especially the following special case: To solve 


ax—1=0 (88) 
we set 
r=1— (89) 
and compute 
xe = Xo + xo(r + 7’). (90) 


We turn now to our example which consists in expanding the reciprocal of the 


Bessel function 
(91) 


J o(z) 


yt (2-4)? (2-4-6)? 


= 
4 
A 
| 
a 
— 
| 
| 
| 
= 
| 
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into a Laurent series between the first two positive roots of this function, which are 
approximately £:=2.4, &=5.5. In order to avoid even exponents we consider 


(92) 
whose reciprocal is to be expanded in Laurent series between its zeros 
t= 5.76 and = 30.25. 


Let us notice that 13 is near the geometric mean of these numbers. In order to realize 
the condition (26) we replace in (92) z by 13z, also changing the sign of the function 
for formal reasons. Thus let 
a(z) = — Jo(4/132) = >> (93) 
n=0 
be the entire function whose reciprocal 
1 
waz” (94) 


Jo(v/ 132) 


we are to expand in a Laurent series convergent on and near the unit circle | z| =1. 
Below are the 10-place values of the coefficients a, of (93) as computed by 


a, = 


| 

| a | | is | A. | l, A, | B, | 
0 | —1.00000 00000 0 | 2.549 122| .000 000] .601 975 
1 | 3.25000 00000 1| 2.262 721 | —.257° 032] 1.063 727| —.361 583 
2 | —2.64062 50000 2| 1.655 481 | —.395 409| .489 993 | —.144 072 
3|  .95355 90278 1.081 379| —.427 381] .219 762] —.062 309 
—.19369 16775 4| .661 333| —.405 462| .097 —.028 189 
s| 02517 99181 5| .379 302 | —.361 625] .042 831] —.012 990 
6| —.00227 31870 6} .193 s00| —.309 968| .018 815 | —.006 018 
7|  .00015 07726 7| .070 603 | —.256 308| .008 261 | —.002 786 
3 | —.00000 76564 |—.011 124| —.202 .003 630| —.001 285 
9/  .00000 03072 9| —.065 027| —.150 752| .001 604 | —.000 588 
10} —.00000 00100 10 | —.099 093 | —.099 722| .000 727| —.000 260 
,00000 00003 11} —.117 947| —.049 615 | .000 —.000 099 
= 39229 24951 12| —.123 985} .000 .000 136 


From these values, rounded to 6 places, we computed to 6 places the values of a(z,) 
at the 24th roots of unity 


= cos (15u)® + 7 sin (15,y)°, = 0, 1,---, 12), 
and from these the values of the reciprocal 


w(Z,) 1/a(z,) R, + (u 0, 1, 12), 
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which are tabulated above. The coefficients A, and B, of the interpolating cosine and 
sine polynomials (75), (76) where then found by the 24-ordinate scheme. They are 
tabulated above. 

From these values we computed the coefficients w” of the approximation 


12 (0) 
wo(z) = wn 2 
—12 
according to (78) by the formulae 
Ao, 
ws = + 
wlan 3(An B,), 
=o = 4d, (m= 11). 


These values rounded to 5 places are in the first column of the following table which 
contains the complete computation according to the relations (89) and (90). The last 
column headed w= +w(r+r?) gives the 9-place values of the coefficients w, of 
(94). 

Remarks. 1. The basic numerical process in this computation is obviously the con- 
volution of sequences. Thus the second column aw is obtained by the convolution 
of the column a@ with the column w™. According to the formula (30) this is done very 
simply rewriting the column a, say, in reverse order, then matching it with the column 
w such that the zero term of one column corresponds to the mth term of the other. 
The accumulated products of matching elements gives the mth term of the product 
column aw), This operation is very familiar from the process of smoothing by means 
of a linear compound formula. 

2. The operation of convolution of sequences implies an important check by means 
of their sums, for it is clear that the sum of the product column should equal the 
product of the sums of the factor sequence, except for the accumulated rounding 
error. At the very bottom of each column we wrote the actual sum of the sequence 
in that column. Directly below it we wrote (in parentheses) the value of this sum in 
terms of the sums of the columns which enter into its composition. 

3. The column of final residuals I—aw was also computed (values not recorded 
here) and its terms were found to be so small that a further repetition of the process, 
with our 10-place values of the a@,, would not alter our 9-place values of the w,. As 
final checks we found by (93) 


a(1)—! = 2.54911 8356, a(—1)-! = — .12398 5065, 
a(i)-! = .19349 9936 — .30996 7383 i. 


The corresponding values of w(z), computed by (94), were found to be 


2.54911 8355, w(—1) = — .12398 5067, 
.19349 9940 — .30996 7383 i. 


w(1) 
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n w aw) r=I—aw r? 
—29 
—28 
—27 
—26 
—25 
—24 49 
—23 4 
—22 —102 
—21 64 
—20 18 
—19 — 64 
—18 85 
—17 — 65 
—16 9 
—15 49 
—14 — 80 
—13 62 
—12 .00007 — .00007 00000 .00007 00000 —- 7 
—i1 .00023 — .00000 25000 .00000 25000 — 68 
—10 .00049 .00007 26562 — .00007 26562 97 
— 9 .00110 — .00004 80946 .00004 80946 — 42 
— 8 .00246 .00002 68539 — .00002 68539 — 35 
— 7 .00552 — .00000 52301 .00000 52301 61 
— 6 .01242 — .00001 62992 .00001 62992 — 27 
— § .02791 .00002 32702 — .00002 32702 — 17 
—4 .06274 — .00001 52800 .00001 52800 24 
— 3 . 14104 — .00000 13045 .00000 13045 — 14 
— 2 .31703 .00001 89326 — .00001 89326 7 
—1 .71266 — .00002 53421 .00002 53421 — 13 
0 .60198 1.00002 07580 — .00002 07580 — .00000 00056 
1 .35107 .00000 39064 — .00000 39064 — 70 
Z .17296 — .00001 57627 -00001 57627 237 
3 .07873 .00000 43626 — .00000 43626 — 89 
4 .03455 .00001 22969 — .00001 22969 —116 
5 .01492 — .00001 61378 .00001 61378 174 
6 .00640 .00000 57617 — .00000 57617 —145 
7 .00274 .00000 04662 — .00000 04662 67 
8 .00117 .00000 78363 — .00000 78363 30 
9 .00051 — .00001 53488 .00001 53488 — 91 
10 .00023 .00001 88876 — .00001 88876 91 
11 .00013 — .00001 23781 .00001 23781 — 37 
12 .00007 .00006 13042 — .00006 13042 — 29 
13 .00002 88471 — .00002 88471 81 
14 — .00009 48797 .00009 48797 — 75 
15 .00004 63585 — .00004 63585 3 
16 — .00001 07391 .00001 07391 53 
17 .00000 14982 — .00000 14982 — 39 
18 — .00000 01411 .00000 01411 — 13 
19 .00000 00096 — .00000 00096 37 
20 — .00000 00005 .00000 00005 — il 
21 — il 
22 0 
23 32 
24 — 21 
25 80 
26 —124 
27 5 
28 103 
29 — 92 
30 43 
31 — 13 
32 3 
> 2.54913 1.00000 45679 — .00000 45679 .00000 00002 
(1.00000 45680) (.00000 00000) 


| 
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n r+r? w(r+r?) =o (7 +r?) 
1 
- 2 
—27 5+ .00000 0001 
—26 12 00000 0001 
—25 26 00000 0003 
—24 49 58 .00000 0006 
4 130 00000 0013 
—102 292 00000 0029 
= 64 657 .00000 0066 
—20 18 1476 .00000 0148 
—19 — 64 3318 .00000 0332 
—18 85 7459 .00000 0746 
i — 65 16767 .00000 1677 
—16 9 37691 .00000 3769 
49 84725- .00000 8472 
—14 — 1 90454 .00001 9045 
62 4 28120 00004 2812 
3 6 99993 2 62369 .00009 6237 
| 24932 —1 36696 .00021 6330 
—10 —7 26465 — 37118 00048 6288 
4 80904 — 68748 00109 3125 
~ § —2 68574 — 27684 00245 7232 
=» 7 52362 36008 .00552 3601 
— 1 62965 — 35203 .01241 6480 
= § —2 32719 9539 .02791 0954 
~§ 1 52824 9186 .06274 0919 
» J 13031 — 49463 .14103 5054 
—1 89319 21184 31703 2118 
= 2 53408 — 48016 .71265 5198 
0 — .00002 07636 — .00000 53032 60197 4697 
1 — 39134 23489 .35107 2349 
2 1 57864 7594 17296 0759 
3 — 43715 — 34341 07872 6566 
4 —1 23085 22715- 03455 2271 
5 1 61552 4594 01492 0459 
6 — 57762 — 20946 .00639 7905 
7 — 4595 — 47102 .00273 5290 
8 — 78333 — 20304 00116 7970 
9 1 53397 —1 15311 00049 8469 
10 —1 88785 —1 73068 00021 2693 
11 1 23744 —3 92531 00009 0747 
12 —6 13071 —3 12836 00003 8716 
13 —2 88390 1 65178 00001 6518 
14 9 48722 70470 00000 7047 
15 —4 63582 30065- 00000 3006 
16 1 07444 12827 .00000 1283 
17 — 15021 5472 .00000 0547 
18 1398 2335- .00000 0233 
19 « 996 00000 0100 
20 425+ .00000 0043 
21 a 2 182 .00000 0018 
22 0 78 .00000 0008 
23 32 33 00000 0003 
24 14 00000 0001 
25 80 6 .00000 0001 
26 — 12% 2 
1 
28 103 1 
29 — 92 
30 43 
31 = 
32 3 
— .00000 45677 — .00001 16443 2.54911 8355 
(—.00001 14637) 
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THE PROPAGATION OF WAVES IN ORTHOTROPIC MEDIA* 


BY 
G. F. CARRIER 
Harvard University 


1. Introduction. The present article is an extension of a previous paper dealing 
with the elasticity problems of orthotropic media.! Here, the displacement potentials 
which define the dynamic phenomena in such media are discussed. 

2. The dynamic problem. Hooke’s law for an orthotropic medium may be written 
in the form 

Or = + + 


(1) 


wherein we use the conventional notation for the stresses and the strains. If we limit 
ourselves to a consideration of those materials which are isotropic in the y, z plane? 
and for which 


bes = b55 = (b11b22 — bi2)/(b11 + bee + 2b12), (1a) 


the number of independent elastic constants is reduced to four and the dynamic 
elasticity problem may be easily treated.* We utilize the familiar equilibrium equa- 


tions 

+ aX = (2) 
Ox oy Oz ot? 


and define the displacements in terms of potentials as 


Ou 
Oy Ox 


When Eggs. (1), (2), and (3) are combined, we obtain 


oy Ox? Ox? dy? 02? or? dz? 

- +2 -0, (4c) 
dz Ox? dy? Ox? dy? 


* Received August 18, 1945. 

1G. F. Carrier, The thermal stress and body force problems of the infinite orthotropic solid, Quart. Appl. 
Math. 2, 31-36 (1944). 

2 The isotropy implies that be =b33, bi2=bi3, bag = —b23)/2. 

3 These conditions are imposed in order that the roots of Eq. (8) will appear in a useful form. They 
include isotropic media as a special case. 
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where a, b, a, y, and 4, are given respectively by bi/p, bee/p, bee/p, (b22 — /2p, 
(a—a)(b—a), and b—y. 

If we now consider the homogeneous equations (i.e., vanishing body forces) and 
require the ¢; to vanish at infinity, we may begin the integration by removing the 
leftmost derivative of each equation, multiplying each term in the remaining forms 
by exp[—i(xt+yn+<f) ], and integrating the equations so obtained over the infinite 
region. We find® 


+ a(n? + + =| Vi + + = 0 (S) 
and two similar equations. Here 
ff $;(x, y, 2, t) exp [— i(xt + yn + 2}) ]dxdyds. (6) 


We now have three ordinary linear differential equations (in ¢) the solutions to which 
are of the form 
¥; = A,(E, §) cos wt. (7) 
In order that these solutions be non-trivial, the determinant of coefficients of Eqs.(5) 
wherein 02/02 has been replaced by —w* must vanish, that is 


ag? + a(n? + — w’, Bn’, Bs? 
Be?, ak? + bn? + yf? — = 0. (8) 
Be, dn’, at? + yn? + bg? — w? 
The three roots of this equation are easily shown to be 
w = af + (9) 
= at +¢), (10) 
w= af ta (11) 
Corresponding to these roots, the A;; must respectively obey the relations® 
Ay:Ay2:A13 = @ — a:B:8B, (9a) 
Ao:A22:A23 = 0:BF?: — Bn’, (10a) 
A31:A32:A33 = (b — a)(n? + §*): — BE: — BE. (11a) 


Because of Eq. (6) it is evident that 


1 
f f tt) exp + yn + 2f) 


k=1,2,3 8x* 


f f f £) Cos wnt exp |i(xt + yn + (12) 


‘ This condition may be justified by noting that functions describing a time dependent phenomenon 
which originated at fy in a defined region of the medium must vanish at infinity for all time ¢. 

5 The procedure thus far has been identically that of the reference in footnote 1. 

6 Here, Ax; is the A; of Eq. (7) corresponding to wx. 


¥ 
— 

I 

ag 
= 
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The solutions of interest are those for which the A,; are explicit functions of w; only. 
We consider first the solution associated with w and find it convenient to write A as 


= (a — a) rary exp 


Here, w: is defined as =ia'/2£+ and By 
is an arbitrary function of n=|n|. 
If 7 is the angle between rm and w,, and r is the polar angle about wi, the value of 


Ai becomes 


Ay = (a — fff Bi,(r1) exp cos sin ududvdr, 


ll 


0 


— (13) 


Using the same notation’ and procedure, we find that @1 becomes® 


ou ff f cos wt exp cos sin 


[(a a)/en] Cy; (w) { sin t) sin w,(r1 + t) \ du 


(29°71) [Bu(ri t) By(ni + t)|. (14) 


Thus we find that one of the possible propagation phenomena is described by a mo- 
tion wherein the wave fronts are the ellipsoids r,=const., and whose amplitude at- 
tenuates like 1/r:1. Equation (9a) indicates that the other displacement potentials 
associated with w, differ from @1 by a constant factor. In fact we could now write the 
displacements for this motion as the gradient of a single potential in a distorted co- 
ordinate system. This fact will be useful later. 

The determination of the motion associated with w, requires only a slight variation 
on the foregcing procedure. As required by Eq. (9b), we define A22= —{?C(w2)/we, 
Ao3=7?C(we)/w:, where C(w:) is derived as before from an arbitrary function B(re) 
and w, and fr are defined in the same manner as were w and rm. The expression for 
becomes 


1 2 — £°C(we) 
oo = os wot exp [i(xt + yn + 2f) 
873 We 


1 C(we) 
= — — fff cos wot exp [i(xt + yn + 2) |dtdndt 


82? dz? We 
1 Biret+t 
| (re (re (15) 
2x* 


Similarly, 


7 The angle v should now be measured about ry. 
8 The final step is due to the inverse relationship implied by Eq. (13). 
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(16) 


The displacement ve associated with this solution can now be written in the form 


1 
V2. = — curl {i 


2x? 


B(re t) 


(17) 


The fact that the equations governing our problem are linear implies that the 
0?/dydz of the foregoing expression is superfluous,® that is, we may write 


1 


2a 


_ — 2b) 
— curl 


(18) 


The displacements corresponding to w; arise in an entirely analogous manner and are 


given by 
i(b — a), i8, kp 
0 0 
ax ay’ as (19) 
0, P(rs — 2) Q(rs — 2) 
r3 rs 


where it is required that 0(P/r3)/d0y+0(Q/rs)/0z=0. The first of the three foregoing 
solutions corresponds mathematically to the potential solutions found for the iso- 
tropic solid and the latter two to the rotational motions; together they suggest a 
way of factoring the equations for the displacements. We write the body forces in 
the following manner, choosing the coefficients of the various derivatives in the opera- 
tors to correspond to those found in the foregoing results; namely 


(X, Y,Z) = grad* (x, y, 2, 2) + curl* M(x, y, 3, 4), (20) 
where 
grad* = i(a — a)d/dx + jBd/dy + kBd/dz 
and 
i(d a), kp 
0 0 
curl* M = — |. (21) 
Ox oy 02 
M,;, M,, M, 


We must require that 0M,/dy+0M,/dz=0. We now assume the displacement in the 


form 
v=grad* ¢ + curl* G, (22) 


substitute into equations (1), (2), and (3), and obtain a set of equations which are 
separable into the following: 


® Actually, all derivatives of B/r2: constitute solutions. 


| 
7 
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a : = — @, (23a) 
Ox? Oy? az? ot? 
0G, dG, 
4 )- (23b) 
Ox? dy? dz? ot? 
| (— 4 0? 4 (M M.) (23 ) 
ax? ay? as? / 


The relation 0G,/dy+0G,/dz=0 will be automatically satisfied. 

Although the above equations may be transformed into the familiar Poisson form 
by trivial transformations, it is interesting to extend the foregoing procedures to ob- 
tain a formal method for the solution of, say, Eq. (23a). The steps leading to Eq. (5) 
transform Eq. (23a) into 


E + ag? + b(n? + T(é, t), (24) 
where y is defined as before and 
T = fff (p, g, s, t) exp [— i(pt + qn + st) |dpdgds. (24a) 
Conventional operational procedures then give a particular integral of Eq. (24) as 
1 t 
y= —f T(é, n, £, sin — a)da (25) 
and ¢ becomes (according to Eq. (12)), 


-exp [i{(x — p)h§+ (y — + da | dpdqds (26) 


1 ed t | 
8x3 


where!” 


| 
I(a) = f f f cos — a) exp [iw -Ri| sin pdudvdw 
0 0 0 


coswi(t—a) . 
ba!/? 0 
Qn sin w;(R; + ¢ — a) + sin — t+ 
f ths (27) 
ba'/?R, 0 @1 


The foregoing integral is a step function which (since 0 Sa St, and Ri >0) has a single 


10 Again we use the type of coordinate transformation which led to Eq. (12). 


ig: 
| 
| 
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step of magnitude at a=t—R;. R; is defined by Ri =(x—p)?/a+(y—q)2/b 
+(z—s)*/b. The evaluation of the Stieltjes integral of Eq. (26) now yields 


1 = t — R 
f f f » apdeds (28) 
Ri 


and we have the familiar retarded potential. The expressions for the components of G 
can be obtained in analogous fashion. 

3. More general media. It is quite evident that one might start with the general 
linear law relating the stresses and strains in an aerolotropic material and by the 
same procedures arrive at three equations analogous to Eqs. (4), using either the dis- 
placements themselves or the potentials defined in the foregoing. One would then ar- 
rive at a determinantal equation of the same form as Eq. (8). The roots could be 
found and Eq. (12) would be valid. However, the w; of this general problem would 
not appear in the concise polynomial form found in the foregoing considerations. In 
fact, the integrand of Eq. (12) becomes sufficiently complex in this general case that 
it does not seem worth while to present in more detail the procedure outlined here. 


aa 
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REFLECTION IN A CORNER FORMED BY 
THREE PLANE MIRRORS* 


BY 
J. L. SYNGE 
The Ohio State Universityt 


1. Introduction. If a plane mirror is attached to the base of a projectile and a 
parallel beam of light projected on it, the direction of the reflected beam at any in- 
stant will give us information about the angular position of the projectile at that in- 
stant. It will not, however, indicate the angular position completely, because a 
rotation of the projectile about the normal to the mirror leaves the direction of 
the reflected ray unaltered. 

This difficulty may be overcome by using more than one mirror, and the possibil- 
ity of using a reflecting corner formed by three plane mirrors suggests itself. If the 
three mirrors are mutually perpendicular, the direction of the reflected beam gives no 
indication of the angular position of the projectile, because such a corner reverses 
the direction of any parallel beam falling on it. But if the angles of the corner are not 
right angles, this is no longer the case; in general, there will be six reflected beams, and 
their directions will determine the angular position of the projectile completely. With 
three parameters at our disposal (the three angles of the corner), we can secure a 
variety of different effects. The purpose of the present paper is to make a systematic 
study of the optical behavior of all corners formed by three plane mirrors. 

The method used is based on the fact that the transformation of ray-directions 
due to reflection in a plane mirror is equivalent to a rigid-body rotation about the 
mirror-normal through two right angles (i.e. a half-turn), combined with a reversal of 
sense. Consequently, three successive reflections in three plane mirrors produce a 
transformation equivalent to three half-turns, combined with a reversal of sense. But, 
by Euler’s theorem, three successive half-turns are equivalent to a single rotation 
(not, in general, itself a half-turn). Thus the transformation due to reflection in a 
corner formed by three plane mirrors may be described by giving the axis of the single 
equivalent rotation (called the optic axis in the present connection), and the angle of 
the rotation. 

It is found that, when different orders of reflection in the three mirrors are taken 
into consideration, there are in general three optic axes and a unique angle of rotation. 
The rotation occurs in both positive and negative senses, so that in general there are 
six reflected rays resulting from a given incident direction. This is, of course, to be 
expected, since we can form six permutations of three mirrors. 

It is useful to represent directions by points on the surface of a unit sphere. There 
are then two fundamental spherical triangles. One has for vertices the normals to the 
three mirrors, and the other has for vertices three directed optic axes. Actually there 
are two spherical triangles formed by (undirected) optic axes, but one is the reflection 
of the other in the center of the sphere, and so we pick out one for definiteness. The 


* Received Feb. 1, 1946. 
+ Part of this work was done by the author at the Ballistic Research Laboratory, Aberdeen Proving 


Ground. 
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vertices of the triangle formed by the mirror-normals lie at the middle points of the 
sides of the triangle formed by the directed optic axes, and the angle of the single 
equivalent rotation is the defect from four right angles of the sum of the angles of the 
triangle formed by the directed optic axes. 

Explicit formulae are given for the construction of the optic axes (3.10) and the 
angle of the single equivalent rotation (4.15). 

The fact that a ray, to be reflected, must strike the front of a mirror, and not the 
back, introduces awkward conditions. These conditions are removed in the mathe- 
matical theory by supposing that the mirrors are planes which reflect from either side. 
Further, in investigating the effect of a second reflection, we may find that after re- 
flection in the first mirror the course of the ray does not bring it into incidence with 
the second mirror. In such cases we shall disregard the position of the ray, and apply 
to its direction the transformation corresponding to reflection in the second mirror. 
These artificialities (from the practical standpoint) are introduced to avoid encumber- 
ing the mathematical theory with conditions which have no bearing on the funda- 
mental transformation problem. Once the general theory has been set up, the condi- 
tions mentioned above may be looked into in any particular case. To facilitate this, 
we shall continue to call one side of each mirror the front. 

2. Equivalence of reflections and rigid body rotations. Let N be a unit vector 
normal to a plane mirror, drawn out from the front (Fig. 1). Let I be a unit vector 
along an incident ray, and I’ a unit vector along the reflected ray. From a point O 
let us drawn the unit vectors N, I, I’, —I, —I’ (Fig. 2). 


It is clear from the law of reflection that I’ is obtained from —I by a half-turn 
about N; equivalently, —I’ is obtained from I by a half-turn about N. 

The transformations may also be shown on the surface of the unit sphere with 
center O; the unit vectors are now represented by points on the surface of the unit 
sphere (Fig. 3). The law of reflection may be stated as follows. Join —I to N by a great 
circle, and produce it on to make the arc (N, I’) equal to the arc (—I, N). Equiva- 
lently, join I to N by a great circle, and produce it on to make the arc (N, —I’) equal 
to the arc (I, N). 

In order that the incident ray may strike the front of the mirror, the arc (—I, N) 
must be less than 37; but in view of the remarks made in Section 1, this restriction 


will not be imposed. 
Consider now a corner formed by three plane mirrors with unit normals Ni, No, 


— 
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= 
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Fic. 1 Fic, 2 _- 
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N;, drawn out from the fronts. These are represented on the unit sphere in Fig. 4. 
Consider successive reflections of an incident ray I in the order N,, Ne, N3. We start 
by marking —I on the unit sphere. We construct the first reflected ray I’ by drawing 


Fic. 3 


a great circular arc from —I through N,, making N;, the point of bisection. If the sec- 

ond reflected ray is denoted by I’’, we construct —I’’ by drawing a great circular arc 
from I’ through Ne, making N, the point of 
bisection. We carry on from —I” similarly 
through N; to form the final ray I’”’. 

We might also have started with I, in- 
stead of —I. The same rules of construction 
would have led us to —I’”’. 

If the reflections occur in a different or- 
der, we change the order of the points 
N,, Ne, N; in the construction. In this way 
we get, in general, six different final direc- 
tions I’’’ corresponding to a single incident 
direction I. 

The transformations described above are 
equivalent to half-turns about diameters 
of a sphere, namely the diameters defined 
by N,, No, Ns, We know by Euler’s theo- 

Fic, 4 rem that any succession of rotations about 

diameters of a sphere is equivalent either to 

no displacement at all or to a rotation through an angle less than 27 about a uniquely 

determined diameter. The former alternative means that I’” coincides with —I, no 

matter how I is chosen. On the other hand, if there is a unique axis of equivalent 

rotation, then rays incident along the axis of that rotation (and such rays alone) will 

undergo reversal as a result of the triple reflection. Every other ray will undergo a 
change of direction determined by application of the equivalent rotation. 

We shall call the axis of the equivalent rotation the undirected optic axis of the 
reflecting corner for the order of reflections assigned. This definition of optic axis 
would be adequate if we were content to have the term denote a diameter, without 
sense of direction. It is, however, desirable to understand by optic axis one of the two 
unit vectors lying on the axis of the equivalent rotation. Accordingly, we shall pro- 
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ceed in the next section to define the directed optic axis. For the present, let us sum 
up our results as follows: 


THEOREM I. For three successive reflections in a given order in three plane mirrors, 
either every incident ray is reversed in direction, 
or there exists a unique undirected optic axis 
such that 

(a) a ray incident along the undirected op- 
tic axis ts reversed in direction; 

(b) the directions of reflected rays are ob- 
tained from the directions of incidents rays re- 
versed by a rigid body rotation through an 
angle less than 2m about the undirected optic 
axis. 


3. Determination of the directed optic 
axes. Consider the following problem in 
spherical geometry: Given three points N,, 
N2, N; on a unit sphere, to construct a spheri- 
cal triangle A;, As, As, such that Ne, Nz; 


are the middle points of its sides (Fig. 5). Fic. 5 
In vector notation, we have 
N, = N,-As;, N.-A; = N2-Ai, = N;-A:, (3.1) 
and also 
= A: + As, L:N, = A; + Ai, L;N; = Ai + A2, (3.2) 
where the L’s are unknown scalars. Our problem is to find the A’s to satisfy (3.1) 
and (3.2). 


Solving (3.2) for the A’s, we get 
2A, IN, + LN, + LNs, 


2A,= — L.N.+ (3.3) 
2A; = + — 
Let us define Mi, M2, M3 by 
M,=N.-N; M:=N3-Ni, M;=N,-Nz, (3.4) 


these being of course the cosines of the angles between the N vectors. Taking the 
scalar products of the first two of (3.3) and N3, we get 


2A, = — 11M2+ 12M, + Ls, 


(3.5) 
Hence, by the last of (3.1), 
L,M2 — L2M, = 0. (3.6) 
Similarly 
— = 0, L3;3M, — 11M; = 0. (3.7) 
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If all the M’s vanish (i.e. if the N vectors are mutually perpendicular), we may 
choose the L’s arbitrarily, save for the condition that the sum of their squares shall 
equal 4. For it is easy to verify that the A’s as given by (3.3) will then be unit vectors, 
and the conditions (3.1) and (3.2) will be satisfied. This means that we can take any 
point A; on the unit sphere and construct the triangle by passing successively through 
the N points, making each a point of bisection. The triangle necessarily closes. 

Now suppose that there is at least one M different from zero. From (3.6) and 
(3.7) it follows that a number & exists such that 


= 2kMi, De = 2kM2, L3 = 2kM3. (3.8) 


We proceed to determine k. Since the N’s and the A’s are unit vectors, if follows from 
(3.3) and (3.8) that 


k M; + M; — 2M,MoM;) = 1. (3.9) 


The quantity in the parentheses is positive definite, since no M can exceed unity in 
absolute value. We have then a choice between two real values of k, one negative and 
the other positive. If either of them be chosen, and the L’s obtained from (3.8), and 
then the A’s from (3.3), the conditions (3.1) and (3.2) are satisfied, and so the prob- 
lem is solved. The problem then admits of two (and only two) solutions; the two 
triangles are the diametrical opposites of one another. 

To avoid confusion, let us pick out one of these two solutions for application to 
the optical problem. Let us decide to take the positive value for k. To sum up, the 
triangle required is given by 


A; = k(— + MWN;), 
A, = k( M,N, — + MSN3), (3.10) 


A; = k( M,N, + M:N, — M3N3), 
where the M’s are given by (3.4) and 
k = (Mi+ Mi+ M3 —2MiMiM;) (3.11) 


the positive value being understood. 
We have also 
N,-A, = N,-A; = 
N2-A; = N2-A; = kM2, (3.12) 


N;-A; = N;-A, = 
and 


2kM\N, = = A;+ Aj, 2kM SN; = A, + Az. (3.13) 


We have now to show the connection of this problem in spherical geometry with 
our optical problem. We shall begin by proving that the diameter through Ae is the 
undirected optic axis for reflections in the order Ni, Ne, Ns. 

We take an incident ray with I= +A». The first reflected ray is then given by 
I’= + A;. The second reflected ray is given by I’’= + Aj. The third or final reflected 
ray is given by I’’ = + Ag. Thus I’ = —I, which proves that the diameter through Az 
is the undirected optic axis for reflections in the order Ni, Nz, Nz. In just the same 
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way we may show that the diameter through A, is the undirected optic axis for re- 
flections in the order N3, Ne, Nu, also. 

Further, if we start from A; or Aj, we can prove in the same way that the diameter 
through A; is the undirected optic axis for reflections in either of the orders Ne, Ns, Ni 
or N,, Ns, No, while the diameter through A, is the undirected optic axis for reflections 
in either of the orders N3, Ni, Ne, or Ne, Ni, Ns. 

We are now in a position to define the directed optic axes by selecting senses on 
the undirected optic axes. We shall do this by imposing the condition that the directed 
optic axes are given by (3.10) with k positive. 

Let us sum up as follows: 


THEOREM II. For any reflecting corner composed of three mirrors which are not all 
mutually perpendicular, the directed optic axes are the three unit vectors A,, As, A3, given 
by (3.10) with positive k. The mirror-normals N,, Ne, Ns meet the unit sphere at the mid- 
dle points of the sides of the spherical triangle formed by Ay, As, As, Ni being on the side 
AcA;, and so on. The directed optic optic axes correspond to the following orders of re- 


flection: 


Directed optic axis Order of reflection 
Ay N;N,N2 or N.NiN; 
As N,N,N; or 
As N.N;N;, or N.N;N: 


We see from (3.10) that the directed optic axes are easily constructed in space by 
adding and subtracting the vectors k MiNi, kM2Ne, kM;N;. This construction is shown 
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in Fig. 6. It is interesting to note that six of the eight vertices of the parallelepiped 
are occupied by the vectors + Ai, + As, + As. 

4. The angle of the equivalent rigid-body rotation. Fig. 7 shows the representa- 
tions on the unit sphere of the three mirror-normals N,, Ne, N; and the three directed 
optic axes Ai, As, A3. We know that three successive half-turns about Ni, Ne, Ns in 
order are equivalent to a rotation about A». We proceed to find the magnitude of this 


rotation. 


As 


A, 


Fic. 7 


We first construct the point B by joining Nz to Ni, and producing the great circle 
so as to make the arc (N,, B) equal to the arc (Ne, N,). We then investigate the dis- 
placement of B resulting from the three successive half-turns. It is clear that the final 
position C is obtained by joining Nz to N3;, and producing the great circle so as to make 
the arc (N;, C) equal to the arc (Ne, N;). The effect of the three successive half-turns 
is to rotate the arc (As, B) into the arc (As, C). Since certain spherical triangles are 
obviously congruent, it is at once seen that the angle of rotation (taken less than 7) 
is equal to the defect from four right angles of the sum of the angles of the spherical 
triangle A, A2A;. If the half-turns were carried out in the reverse order, then C would 
go to B; the angle of rotation would be the same in magnitude, but would have the 
opposite sense. 

If the other senses of application of the three half-turns are considered, we get 
the same magnitude for the angle of the single equivalent rotation in all cases. 

We shall now obtain a formula for the angle of the single equivalent rotation in 
terms of the angles between the normals of the three mirrors. 

If an incident ray I is reflected in a mirror with normal N, the reflected ray I’ is 
given by 
I’ =I — 2N(N-]). (4.1) 
Let us start with an incident ray I such that —I=B, where B is as in Fig. 7, and let us 
follow this ray through successive reflections in the mirrors with normals N,, Ne, Ns. 
The first reflected ray is I’=Nbe, and so by (4.1) 


I =I — 2N,(N,-I’) = N, — 2M;N,. (4.2) 
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The second reflected ray is I’’= —Ng¢. The third and final reflected ray is 
2N;(N;-I’’) a +. 2N;(N;-N2) + 2MN3;, (4.3) 


the M’s being as in (3.4). 
We note that (3.10) give 
= = kM». (4.4) 
As already seen, the angle of rotation 6423 is the angle between the arcs (Az, B) 
and (As, C), as shown in Fig. 7. This is the same as the angle between the vectors 
(A: XB) and (A: XC). If the rotation is considered positive when it is a right-handed 
rotation about A», through an angle less than 7, then the angle satisfies the equation 


(Ae X B) X (Az X C) 


sin = lA, XB||A,XC| (4.5) 
where B= —I and C=I’”. (The angle 623 shown in Fig. 7 is negative in the defined 
“— evaluate the right-hand side of (4.5), we note that by (4.4) we have 

| Ae XB||A,xXC| =1— (4.6) 
Also, identically, 
(A> X B) X (As XC) = X C)]. (4.7) 


By (4.2) and (4.3), 
BxC 


—IXI” = (—N: + X (— + 
— 2M,(N2 X N;) — 2M;3(Ni X — 4M\M;(N3s X Ni), (4.8) 


and so, by (3.10), 
A,:(B XC) = k(M,N, — + MSN;)-(B XC) 


= kP(— 2M; + 4M,M2M; — 2M3), (4.9) 
where P is defined by 


P= N,-(Nz x N;). (4.10) 
When the value (3.11) for k is used, (4.9) may be written 
A,-(B XC) = — 2Pk (1 — k (4.11) 


On substitution of this expression in (4.7), and then substitution from (4.6) and (4.7) 
in (4.5), we get 
sin 0323 = — 2Pk. (4.12) 
We note that P is the determinant formed from the direction cosines of the three 
normals to the mirrors, and 

1 M; 

P? =| Ms; 1 M,|=1-— (4.13) 
M, M, 1 


Hence 
P =e(1 — k-*)'/?, (4.14) 
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where ¢ is +1 or —1 according as the orientation of the triad Ni, Ne, N3 is positive 


or negative. 
We may sum up as follows: 


THEOREM III. The magnitude (taken less than 7) of the angle of the single equivalent 
rotation is the defect from four right angles of the angle-sum of the spherical triangle on 
the unit sphere whose vertices represent the directed optic axes. If a positive angle corre- 
sponds to a positive (right-handed) rotation about the directed optic axis involved, the 
angle is given both in magnitude and sign by 


sin = — 2ek-"(1 — k-2)1/2, (4.15) 


where k is given in terms of the cosines of the angles between the mirror-normals by (3.11), 
and €is +1 or —1 according as the triad of mirror-normals, in the order of the reflections, 
is positive or negative (1.e. right-handed or left-handed). 

Let us now see how the six reflected rays are to be constructed when the incident 
ray Lis given and the three directed optic axes are known. We mark on the unit sphere 
(Fig. 8) the directed optic axes A;, As, As, and the incident ray reversed, —I. Let us 


Fic. 8 


suppose that the mirror-normals are so numbered that Ni, Nez, N; form a positive 
triad. To obtain the reflected ray resulting from reflections in the order N,, No, Nz 
we use the optic axis As. We draw the arc (Az, —1), and rotate it about A; in the nega- 
tive sense through an angle equal to the defect from four right angles of the angle-sum 
of the spherical triangle A; or equivalently through an angle sin~! 2k—1(1 — 
as given in (4.15). This gives the reflected ray Ij23, the subscripts indicating the order 
of the reflections. To obtain Ij3{, we rotate the arc (Az, —I) about A: in the opposite 
sense through the same angle. 

Using the other optic axes, we obtain similarly the whole set of six reflected rays. 
These are shown in Fig. 9, the great circular arcs being shown as straight lines for 
simplicity. All the marked angles are equal. 
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5. Cases of perpendicular mirrors. Let us consider the case where two of the 
mirrors are perpendicular to one another. Let us take N, perpendicular to Ne, and 
write 

M, ~ 0, M2 ¥ 0, M; = 0. (5.1) 


Then (3.10) give for the directed optic axes 


Ai k(— MN, + MN;), 
A> k( MN, MN,), (5.2) 
A; = k( MN, + M.N,), 


; and from (3.11) we have 


3 
0 
K M,N, 
N3 
M,=0 
Az 
Fic. 11 
k = (Mi+ Mi)”. (5.3) 


The equations (3.12) and (3.13) become 
= N, -A; = N,-A; = N.-A; = kMo, N;:Ai = N; = 0, (5.4) 
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and 
A.+ As = 2kM Ni, A; + A, 2kM A, + Ap = 0. (5.5) 
We note that A; and A: are opposed to one another, and Aj, As, A3, N,, Ne are co- 
planar. The arrangement on the unit sphere is shown in Fig. 10. The construction in 
space for the directed optic axes is shown in Fig. 11. 
Let us now consider the case where one of the mirrors is perpendicular to the two 
others. Let N; be perpendicular to Nz and N;. We write 


M:~0, M:=0, M;=0. (5.6) 
From (3.10) we get for the directed optic axes 
Ai=—kM.N;, A:=kM.Ni:, A; = kM\N,, (5.7) 
where k= | M;| —l; from (3.12) and (3.13) 
=N,-A; = kM,, N2-A; = = N;-Ai = N;-A; = 0, (5.8) 


and 

A.+ = 2kMN,, A; + A, = 0, A,+ A, = 0. (5.9) 

N,=As=A3 N\=A, 
N, $1 
$1 
N Ns 
3 Az=A3 in 
tr 
A, 
M,>0,M=M3=0 M,<0,M,=M3=0 
Fic. 12 


Now A; and A; coincide, and A; is opposed to them. The arrangement on the unit 
sphere is shown in Fig. 12 for the cases M,>0 and M, <0. 


177 


ON GRAEFFE’S METHOD FOR SOLVING 
ALGEBRAIC EQUATIONS* 


BY 


E. BODEWIG 
The Hague 


In the usual descriptions of the methods of solving numerical algebraic equations, 
Graeffe’s method takes a minor place as compared with the methods of Newton, 
Horner, and others. It is not useful, of course, for correcting a single approximate 
value, as the other methods are, but has the advantage that no first approximation 
need be known. A second advantage is that approximations to all roots are obtained 
simultaneously, in contradistinction to the other methods which furnish approxima- 
tions to one root at a time. In spite of this, the computations required by Graeffe’s 
method are not much more laborious than those necessary to obtain an approxima- 
tion to a single root by one of the other methods if allowance is made for the time 
necessary to find the first approximation. Yet this slight increase in labor may be the 
reason that Graeffe’s method is somewhat neglected. Its third and perhaps its main 
advantage is that it also affords a means of finding the complex roots. It is true that 
by certain other methods, such as that of Newton, an approximation to a complex 
root can be improved, but obtaining the first approximation is rather difficult in the 
case of complex roots. A last advantage of Graeffe’s method is that it automatically 
separates roots which are close together, such as »/5/2 and 3/2. It is known that 
Lagrange claimed this same advantage for this method of developing a root into a 
continued fraction, in contradistinction to Newton’s method; Lagrange’s method 
fails, however, in the case of complex roots. 

These advantages are well known, though not sufficiently appreciated in practice. 
But so far as can be seen, it is not known that Graeffe’s method also gives the multiple 
(real or complex) roots, in a manner essentially simpler than is generally pointed out in 
more elaborate descriptions of the procedure. Further, of all the methods it is the only 
one for solving an equation having several pairs of complex roots of the same modulus. 

To derive these and other properties, for example the convergence of the process, 
we discuss the whole method in a somewhat simpler form than is generally used. 

Preliminary remarks. The method consists in deriving from an equation 


+ (1) 
with the roots x1, x2, -- - , X, another equation 
(2) 


having the roots X,=x?, where p is a large number (and for practical reasons a power 
of 2), so that the distinct roots of (2) are widely separated and can thus be easily 
calculated in the following manner. 

Splitting of equation (2). 

First case. All the roots of (2) are positive and simple. 


* Received Sept. 28, 1945. 
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Then, from 


we have 


— As = = X1X0Xs ete., 
since the first members of the sums predominate. Thus 
X, = — Aj/1, X2 = — A2/A1, X3 
that is, equation (2) is split into the » approximate linear equations 
X +A; =0, + A2=0, + A3 = AniX + An = 0. 
Second case. Several of the roots of (2) have the same absolute value. 
(ia) There is one pair of complex roots, say 
X3 = Ret, Xe = Re“, 

so that 
Az = XiX; 
= + = X1X2(X3 +X4) = X1X2-2R eos A 
X1X2X3X4 = X1X2: R’. 


Thus X3, X4 are approximately the roots of the quadratic equation 
A2X? + A3X + Ag = 0, (3) 
so that equation (2) is split into the »—2 linear equations 
X +A, =0, AiX + Ap = 0, + As , AniX +A, = 0 


and the quadratic (3). 
(ib) There are two pairs of complex roots, having the same absolute value R, 


= Reta, X5.6 = Ret, 


Then 
Ag = 
— As = + X4+ X5+ Xo) = X1X2:2R(cos A + cos B) 
Ag = + X5X6) = X1X2-2R2(1 + 2 cos A cos B) 
— As = X\X2(X3X4X5 +--+ + X4X5X6) = X1X2-2RX(cos A + cos B) 
Ag = = X1X2-R*. 
Therefore X3, - - , X¢ satisfy the equation 


AoX4 + + AgX? + AsX + Ao = 0. (4) 
To solve it, let us first compute R from 
= (S) 


—-A,=) Xi= 
Ag = >, XiX; = 
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With this R let 

X = R-Y, where Y = (6) 
Then all the roots of Eq. (4) in Y have the modulus 1, and if Y is a root, then 1/Y 
is also a root, that is, Eq. (4) is “reciprocal” : 


BY*+CY?+ BY +1=0. 
By 
Y+Y"=2cos¢=Z (7) 


it becomes a quadratic equation in Z. 
(ic) There are » distinct pairs of complex roots having the same modulus R. 
Then in a manner similar to (ib) above an equation M of degree 2y is split off. 
The modulus R is obtained from the equation 


R** = quotient of the last coefficient to the first coefficient Az of the equation M. (8) 
To compute the arguments A, B, - - - of the uw pairs of roots, we again set X = RY 


as in (6), and divide by the leading coefficient A2R**. The new equation in Y is again 
reciprocal: 

+1) + 4+ VY) +C(¥*?+ ¥2)+---+ =0. (9) 
By substituting (7) into (9) it can be seen that (9) is an equation in Z of degree yu 
having u real roots of Z<2 which can be found by Graeffe’s method. Complications 


cannot arise since all the roots are distinct. The transformed equation therefore will 


break up into a system of u linear equations. 
(iia) There are three roots with the same modulus R, one positive and the others 


complex, say 
X3 = R, X45 = Rett, 


Then 
= 
— As XyX2(X3 + Xa + Xs + Xo) = Xi1X2-R(1 + 2 cos A) 
Ag = XaX5 + XaXs) = X1X2-R*(1 + 2 cos A) 
— As = X1X2X3XuX5 = X1X2: 


so that X3, X4, Xs are approximately the roots of the cubic equation 


AoX* + A3X* + AqX + As = O. (10) 
The value of R is obtained from 
= — A;/A2, 
or more simply from 
R= — A4/As. (11) 


Thus equation (2) is split into »—3 linear equations and the cubic equation (10). 
The last equation can be broken up into the linear equation (11), that is AsX+A.=0, 


and a quadratic equation. 
(iib) There are » pairs of complex roots and one positive root, all of the modu- 


lus R. 
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An equation M of degree 2u+1 results and from that a linear equation L is again 
split off. It consists of the two middle terms of M, that is 


L: Ayi2X + Ayys = 0 (12) 
with the solution R. 

To find the arguments of the complex roots, we again set X = RY and obtain a 
reciprocal equation in Y of degree 2u+1. It has the solution Y=1, so that by dividing 
it by Y—1 there results a reciprocal equation of degree yw in Z that may be solved by 
Graeffe’s method, yielding u pairs of complex roots. 


(iii) There are multiple roots. 

Let X2 have the multiplicity v. Then equation M, mentioned above in (ic) and 
(iib), will be divisible by (X — =(X —R)’, if Xz is real, and by [(X —X:)(X — X2) |’ 
=(X*—2R cos AX+R?)’, if X2 is complex. The reciprocal equation in Y is then di- 
visible by (Y—1)” or (Y?—2 cos AY+1)’, respectively. 

Note. It is not always possible to eliminate the multiple roots of (2) by eliminating 
at once the multiple roots of (1), for distinct roots of (1) may, by the successive squar- 
ings, give the same roots of (2). 

In summary we can say, if the absolute values of the roots of (2) are partly equal, 
partly different, then Eq. (2) is split up into several approximate equations M, of 
lower degrees. The degree of an M is equal to the number of roots having the same 
modulus R. Thus there are as many equations M as there are distinct moduli. To 
every simple root there corresponds a linear equation. 

Determination of the equations /;. It is well known that by squaring the roots of 
(1) a series of equations Gi, Go, Gs, - - - having the roots x?, xf, xf, --- results. 

The problem now is to decide which equation G first breaks up into equations M,, 
and what these M; are. 

We have seen that if the equation M has m roots with equal moduli R, then the 
absolute member of the normalized M (that is an M whose leading coefficient is 1) 
is equal to +R”. In the following transformed equation it will be equal to +R, 
that is, from a certain equation G;, the absolute member of a (normalized) M is (to 
the required degree of accuracy) squared when passing to the following transformed 
equation. This or a similar relation does not hold for the other coefficients of M, for 
they involve cos A. Since cos A changes to cos 2A in the following equations, the co- 
efficients not only irregularly change their quantity, but often their signs too. 

To find the various M; into which an equation G;, is eventually split up, we must 
therefore seek only those coefficients that are squared when passing to the next equa- 
tion G. That is, we begin with the leading coefficient 1 of all G, and choose the first 
member A; that by the last root-squaring is itself squared. The coefficients from 1 
to A; form the first equation M,. The next equation Mz has the leading coefficient A. 
Since A; is itself squared when going a step further, it is unnecessary to normalize the 
supposed equation M2, but we may choose immediately the first coefficient A ; after A; 
that is squared by the root-squaring. Equation M; now has the coefficients lying 
between A; and A ;. If A; is the next coefficient after A ; that is squared by root-squar- 
ing, M; extends from A; to Ax, and so on. 

It is not necessary that the same equation G, yield all the various M;. The higher 
the degree of an M is, the later the mentioned quality of the extreme coefficients will 
generally appear. However, from the stage where an M is split off, it is no longer neces- 
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sary to keep it during the further calculation. It is sufficient to treat only the rest of 
the G, that remain after cancelling the M;. 

Resolution of an equation M/. Every M of degree m must be solved separately, 
in the following manner. 

(i) Normalize M to M’. 

(ii) Find the modulus R of all roots of M from (8), that is, in the normalized M’, 


from 
R™=—absolute member of M’ if m is even, 


and from the linear equation L (12), if m is odd 

(iii) With this R set X = RY, where Y=e**, into M’ and normalize again to M”’. 
This new equation in Y is reciprocal. 

(iv) If this is possible, divide M” by Y—1 (that is, is Y=1 a root?), and repeat 
this division as often as possible. If this is possible s times, then M has the root R 
of multiplicity s. If m is odd, s is at least 1. : 

(v) Form the quotient Q= M’’/(Y—1)*; this is also a reciprocal polynomial. 

(vi) In Q=0 set Z= Y+Y~. Then Q=0 is transformed into an equation Q’ =0 
in Z of degree (m—s)/2. The equation Q’ =0 has all its roots real. 

(vii) cos 6=Z/2 yields (m—s)/2 values of ¢ and therefore m—s values of X: 
X = Re*‘*, and this together with the root X = R of multiplicity s yields the complete 
system of roots of M. 

Solution of equation (1). To every root X; of M there corresponds one root x; 
of (1). Since X;=27, every X; yields p tentative values of x; from which the right 
root must be chosen. We do not, of course, calculate with the complex values of x;, 
but take only the real component of Eq. (1) that is: 


R* cos + R*a, cos (n — + cos (n — +a, = 0, (13) 


where ¢ is given by 
= (2gx + ¢1)/P, q=0,1,2,---,p-—1. (14) 


After having found the first ¢ satisfying (13), we stop the calculation with this ¢;. 

If equation (1) is not too complicated, that is, if the M’s have a low degree, the 
process can be abbreviated in the usual manner. 

In all other cases the process of finding the complex roots can always be abbrevi- 
ated in the following manner. In the chain of the transformed equations G; we go back 
to the first equation G,,(m<k), where M starts to split off, that is, where the double 
products have no more influence on the first (or second) decimal place of the corner- 
coefficients of M. This equation M is solved to one or two decimal places only. The 
root of Eq. (1) is now—to one or two decimal places—to be selected from a group of 2 
members, instead of from the 2* members of (14). This decreases the work involved. 

In addition, the process is abbreviated by the fact that not all 2* equations of (13) 
have to be computed, since the members of the majority of them differ from those 
of the others only in the sign. An example will make this clearer. 

Nevertheless this part of the labor is the most tedious of the whole process if 
the first transformed equation Gm, from which M is first split, has a large m, that is, 
if two or more roots of (1) are close together. It is good, however, to have a method 
that yields these roots at all. 
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In cases where there are only one pair or two pairs of complex roots this whole 
process is superfluous as will be seen in the example. 

Convergence of the process. If the coefficients of the transformed equation G; are 
determined exactly and the roots of G; are calculated by splitting G; into several M, 
then the roots of the M’s are only approximately the roots of G,. If, thereby, a root 
of G, is found with the relative error €,, then the error of the same root in the follow- 
ing equation is 

~ 
as is easily seen. Let us suppose that we are dealing with the largest root X, and there 
are m complex pairs of roots of the same modulus R, then we see that R is calculated 
from the absolute member a of the equation M;. Let now X:=R-10~ be the follow- 
ing root and let the m pairs of complex roots be 


X’ = Reta, Xx" = RetiB,---, X(™ = Rett, 
Then a, being the sum of the combinations of all roots of G; by 2m, is equal to 
a = + A+ cos B+ --- + 
= R?(1 + 2S-10-‘), 


where S is the sum of the cosines. Thus 
2m 
R = Va (1 — (1/m)10-'S), 
that is, since S<m, the relative error e, of R is less than 10+. 
From the following equation G;41, equation M,,4: is set up, and if b is the absolute 
member of Mx4:, then 


b= > where < m. 
Thus 


R’ = R? = (1 — 10-9) 


That is, the relative error of R’ is 
—2i 2 
€k+41 < 10 

Now the roots x1, x2, x3,:-- of (1) are the 2*th roots of the roots X; of G;: 
Thus if € is the relative error of X; of Gx, then 

2k 
m= = (1+ €/2*). 
That is, the relative error of x; is €/2*. But from G41 we obtain, as we have seen 
= = VX1-(1 + 
that is, the relative error of x; is e?/2*+! or the square of that of the equation G,;. Thus, 
we have: 

The relative error of the roots x; of (1) as they are computed from the transformed 
equations G; decreases quadratically with every following equation G,, that is, if the roots 
x; of (1) following from G;, are exact to r decimals, then equation G;+1 will yield them to 2r 
decimals exactly. Roughly, every following equation yields tiwce as many exact decimals of 
the roots x; of (1). 
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It must be taken into account that this property holds only if the roots are already 
sufficiently separated, for instance, if the difference of any two neighboring roots 
of G, is at least equal to 100, or else the approximations above are invalid. 

From this property it follows that Graeffe’s method has its greatest efficiency if it 
is carried out to many decimal places. If a calculating machine is used this does not 
require more computational work than required for fewer decimal places. Now, it is 
true that in this case one must calculate more transformed equations at least if the 
number of decimals are to be fully used. But for this purpose it will be sufficient to 
have one or two equations more. If, for instance, the equation G; yields 5 exact decimals 
of the x; of (1) then the next smaller root of Gs; has a modulus r~10-5R, where R is 
the modulus of the greatest root. In G¢ the relation of the two greatest roots will then 
be 10~"°, that is, only the 10th decimal place of the coefficients of Ge will be influenced. 

From this it follows that—apart from exceptional cases—the same number of trans- 
formed equations will in general be necessary if a certain exactness 1s required. For, 
suppose we have an equation with two roots having the ratio 1.1. Then by 3 or 4 
transformations the ratio will become 3. Thus to have a certain exactness, it will 
be necessary to calculate 3 or 4 equations more than for an original equation with two 
roots having the ratio 3. The example is very unfavorable, for there will be few equa- 
tions having roots of the ratio 1.1. At a ratio 1.5, there are only one or two more 
transformations required. 

Since by raising to powers all roots with distinct moduli will be separated auto- 
matically by a quadratically convergent process, Graeffe’s method is more powerful 
than other ones. 

Influence of rounding off. These considerations of convergence are strictly valid 
only if the coefficients of the transformed equations are exact, but in reality the cal- 
culation is carried out to a fixed number v of decimals. The errors of rounding off are, 
of course, increased by every squaring and multiplication. Now, these errors can be 
estimated by adding the proper inequality to every coefficient of the scheme. But in 
general this tedious supplement will be superfluous, for on the whole the error will 
be annulled by the process of extracting the 2*th roots of the roots of G,. That is, 
if the calculation is carried out to v decimals, the roots of (1) will on the whole be exact to 
v decimals too. 

Moduli or roots lying close together. When the equations of the chain do not soon 
show signs of an approaching splitting up, this will signify that some moduli of even 
roots are close to each other. 

Various procedures have been proposed for accelerating the convergence. But 
they are all unpractical. For they require much more work than does the Graeffe’s 
method when carried on two or three steps further. Add to this that by these de- 
vices the calculation of the chain is interrupted which is very undesirable when at 
the end of the calculation the roots of (1) are to be computed from those of the M;. 

These procedures make no allowance for the quadratic convergence (or diver- 
gence) of Graeffe’s method. For before it becomes obvious that several moduli or roots 
are close to each other several transformations are already effected. By then the roots 
will be separated so far that the greatest difficulties have been overcome, and the 
further calculation will proceed rapidly. It is not advisable, therefore, to disregard 
the entire calculation performed up to this point and instead, to apply Graeffe’s method 
to a transform of the original equation. 
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For instance Encke (Gesammelte mathematische und astronomische Abhand- 
lungen, vol. 1, Berlin, 1888, p. 185), when dealing with a certain equation states 
that “after 6 or 7 operations we have got the conviction that two trinomial factors 
lying close together are existing here.” Then he abandons Graeffe’s method and starts 
on a new Calculation. Yet, if Graeffe’s method is carried 2 or 3 steps further, all roots 
separate automatically. 

Thus the usual procedure of Graeffe will be always the most suttable one. Indeed, one 
of the chief advantages of the method is, that no special devices are required. 

On the other hand, if some knowledge of the position of the roots of f(x) =0 is 
not furnished by the first steps of Graeffe’s method, but by other sources, then this 
knowledge may be used from the beginning to accelerate the separation of the roots 
by transforming the equation first. 

If, for instance, several moduli are close to each other and their absolute value p 
is known approximately and if also several roots are close to each other and their 
values, too, are known approximately, we may proceed as follows. 

We have a group G of roots lying near a circle C with radius p around the origin 
of the Gaussian plane. And in this group G there exist several places u, v, w, - 
where the roots “accumulate” so that.G is divided into the subgroups U, V, W,---. 

Now the slow convergence of G when applying Graeffe’s method arises from the 
fact that the quotient of two moduli of G is lying close to 1. This difficulty will be 
partly overcome by chosing the origin of the coordinates in the neighbourhood of one 
of the points u, v, w, « - - , say u. The circle C will still be a circle, but it does no longer 
have the new origin as its center. Thus the quotients of the moduli of the group G 
have essentially changed. For the relations of the distances of the subgroup U from 
the origin to each other as well as to those of the subgroups V, W,--- differ now 
widely from 1. The same holds for the relations of the distances of the group V to 
the distances of W, - - - . However the relations of the distances of V to each other 
remain nearly the same as they were before, and the same will hold for W to each 
other. 

By this transformation the separation of the group G will, therefore, be acceler- 
ated, that is, the equation of the group G will break up into the equations for the sub- 
groups U, V, W, - - - faster than would have been the case without this transformation, 
and the subgroup U will even be split into its individual elements. The method is 
to be repeated eventually, as far as the subgroups V, W, - - - are concerned. 

The value of this method is largely theoretical because equations with several 
points of accumulation, u, v, w, - - - , near the same circle C do not occur frequently. 

The details of the transformation mentioned above will depend on the nature of 
the roots: 

4. Several roots near the circle C are real. We may suppose that all these real roots 
are positive. Otherwise we form the first transformed equation and get a new equation 
with the assumed property. The convergence will be most rapid if the origin of the 
coordinates is chosen in the neighborhood of the least of these positive roots, say a, 
i.e., if the following transformation is made: 


x=y+a’, where a’ 


ut. The roots near C are all complex, but “simple,” i.e., no two of them are close to- 
gether. In this case, too, the transformation 
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«= y+p’, where p’ ~p, 


will be sufficient as is evident geometrically. 

vit. All roots are complex, but several of them are lying near u=a+1ib (and u’ =a—1b, 
of course). Then it will not be convenient in general to bring u into the origin immedi- 
ately, for this would require a complex transformation. It will be more suitable to do 
this by two real transformations. First we make the transformation 


x= y+a’, where a’ <a, 


so that uw will come near the y-axis. The equation in y will have several roots near +7). 
We transform it by Graeffe’s procedure and obtain an equation in Y having several 
roots near the point (0; —b?). Then by the second transformation 


Y=U-B 


we obtain an equation in U having several roots near the origin. The ratios of the 
moduli will now differ widely from 1. 

This method will be particularly useful when all the roots of the equation are lying 
near the circle C, for instance, for the equations M;. In this special case the method 
may be brought into a more convenient form by applying a procedure of Ostrowski. 
[Recherches sur la méthode de Graeffe et les zéros des polynomes et des séries de 
Laurent, Acta mathematica, 72, 245 (1940) ]. 

In case iii above, i.e., if all roots of f(x)=0 are lying around u=a+1ib and 
u’ =a—ib, the sum ma of all roots will be equal to the coefficient —a,, so that the real 
part of all roots will be approximately @~ —a,/m. After the transformation x =y+a, 
the coefficient of x™~! will then vanish. The same holds for the equation in Y. That is: 

If we know that f(x)=0 has all its roots near two conjugate complex numbers, 
we may apply the transformations of iii without knowing these points, by bringing 
f(x) into the reduced form f(x), transforming f once by Graeffe’s procedure into F(x) 
and bringing the latter into its reduced form F(x). The :0ts of F(x) =0 differ widely, 
and Graeffe’s method will converge then rapidly. 

Exampie. By the following example we skow tk: «fciency of the method in the 
case of roots lying close together or having the same » ¢1i1li. 


825 + 4x4 + 18x”° — 15x? —- 18% — 81 = 0 or normalized: (A) 
x + 0.54 + 2.25x* — 1.875x? — 2.25% — 10.125 = 0. 


In the coefficients of the transformed equatiotis there is always a power of 10 
omitted; its exponent is given in italics on the left of the coefficient, so that, for in- 
stance, the last coefficient of G2 is in reality —104-1.0509---. 

To avoid slips it is advisable to put under each transformed equation G;, the signs 
of the equations having the same roots, but of opposite signs. These signs are alter- 
nately equal or opposite to the signs of Gx. 

As may be seen, it is not until G; that the pace of the coefficients begins to become 
more regular. This late start indicates that the moduli of the roots of (A) lie close 
together. Also from this point of view it is advisable to carry out the calculation to 
many decimals. Equation G; is the first equation where the approaching split is 
perceptible, for in the double products forming the coefficient of x* two zeros appear. 
From G; onwards the process goes rapidly, so that in Gs the coefficients of x, x’, x® 


BS 
Sty 


[Vol. IV, No. 2 


E. BODEWIG 


186 


6 8L860 OSOP‘Z— 


£82 


$ SZLSE £9000 0000‘ 0— 
9 PHLLL LSIIL 6— L BOLL*L 9£982 O£100 0000‘0— 8 9ZOLE OBSE‘Z 
6 SLZ78Z 9100'0— 6 Z8LZ9 61000 0000‘0— 
L=4 LZ6€8 ZS9ZZ BI9E‘ EZ 8 FILET SO9BS 0 £8080 94606 £000‘0 Z ¥SSZE 7998'9 
8 66916 Z1£9'6 9 $4990 — £ 9B8SLZ 19 9 $7696 SEEL‘EI— OF 
‘ + + + 
£ 69019 OlbZE — 9 SE6LZ 8 61979 It T 6 £0687 £6198 = ST 
L £0791 69S7‘0— L L86L8 ILIST ¥000'0 
Z8b6S BOSOE £ 66£6S 8SL0'Z 92 L OOT9L «OZ 8 SSS8E $I O 6PLOS Z 
££ OL076 6LE8Z E9SE‘O— Z 100S8 6SZL6 Z160'0 
s=4 OS160 FLEIT Z9L8°1— 99 96186 SSOL‘S § 90099 PE O8ZLO Z1S66 S697‘ 
L8 SIO£O FIFL6 LE80'Z 42 OF LETO‘S— 12 SL£8€9 ZOELO ST SS6%6 
9 I6Z0L LbS68 861Z‘I— 9 0900S 76679 6676L 90SL8 0069'L—- = OF L 09826 Z 9 SILI 
61 991S‘0— ZO OSSSF 19£66 
$8 YEETL £Z $9S96 EZLET SL 9OZFE LLZ6S Z8ZSO 8600‘F 
06 ZHEL6 O8LZE tI $6 HIZEL 8Z61°9— IT 06 ZZLLO 9610'F 0 91F9L 07165 
+ + + + + 
Z WISE 10198 = 8 T OSLLZ 20896 Z L 6SOI8 F698. L 8£967 16920 6b00'Z 48149 2 
£=4 £ 67908 16£8‘9FT SZELL PELOZ 9 S9TZL ILSZ‘ $Z9S9 166S‘0— 
0 8ZSPE 86000 SISIS 86876 9 1016Z EZ10L 6766'8+ $ 79S10 2 
SZ 9OFI6 6OSO‘I- & 0 OSZIS 0079‘ 2 SLE60 9bZZS 0986'9+ 2 SZ18Z 8866°7- TT 
© SZI8Z S 18962 OSZI8 Z886'Z 
0 0SZ90 68Z1Z 8Z80'T £ 0 SZ790F 16196 Z 
‘ ‘ + + ‘ 
osz9s Iszo‘I- 2 sz 0 0 
Oszt‘or 
sz90'S+ 0 sz 0 $z90'S 0 sz‘o- 
0 00000 00000 $zt0‘t— 0 00000 00000 0. © 00000 00000 0 © 00000 00000 © 00000 00000 0 


= 
| 
| 
| 
| | 
| 
| | 
| 
| 


1946] ON GRAEFFE’S METHOD 187 


are fixed to 16 places since they are no longer influenced by the double products. 
Thus Gs is split into the two equations M: 


M, = X?* — 2.358036915546003 -10°X + 1.390084523771616-10'? = 0 
= 10'*-1.390084523771616X* — 10'*7- 1.287533728649992 X? 
+ 10?!?-1.545884788849941X — 10%*-2.405072447095789 = 0. 


Il 


These two equations must now be solved. Since in Gs the coefficient of X‘ is ap- 
proximately twice as large as the coefficient of X* in G;, this signifies that M, has a 
real double-root, so that M, may be put into the form 


M, = (X — = X? — 2mX + m’? 

as is approximately confirmed. Thus from the coefficient of X we have 
X, = X2 = 1.17901845777300- 10%, 

whereas from the coefficient of X° we have 
X, = = 1.17901845777393 - 10". 


The difference of the two values arises from the rounding off errors we mentioned 
above. To annul them slightly we could take the mid-value of the two values: 


X, = = 1.17901845777346 - 10%. 


Equation M: no longer splits into factors as can be seen from the course of 
its coefficients. Otherwise some of the double products should begin to converge to 
zero. Now, this is only the case with the product 2A4Ao, and there only because of the 
coefficient A,. This follows from the fact that equation M, is already split off and has 
no influence on Mz. Thus Mz has all its roots with the same modulus R. Since its de- 
gree is odd, we have according to p. 179, ii and according to (12): 

R = 10*!*-1.545884788849941 /10'® 1.287533728649892 


= 10-1.200500425311787. 
We now set X = RY in 1/2 and get the reciprocal equation in Y 
M” = BY® — 1.855595246409359Y? + 1.855595246407256Y — B= 0 


or, to make the equation wholly reciprocal, instead of the two middle coefficients we 
put their arithmetic mean, and divide the resulting equation by Y—1: 


Q = M”/(Y — 1) = Y? + 0.2284659663166439Y + 1 = 0. 


Since this equation is no longer divisible by ‘Y—1, we put Z= Y+ Y~ and obtain 


instead of Q: 
Z = 2 cos ¢ = — 0.2284659663166439, 


@ = 2929.71149255575355. 


(The notation “g” refers to the division of the quadrant into 100 parts, instead of 
into 90.) Thus the roots of M2 are 


X3 = R, = Ret‘, 
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The roots of the original equation arise from the X; by extracting the 256th root. Now 


vi. su vl 13 
$/256 = 19.143405439670912305, 


so that the roots x; are to be found among the values 


V3-(cos 2kr/256 + i sin 2kw/256) 
and 
1.5(cos + isin dx), where = + ¢)/256, k =1,2,---, 256. 
There are two roots of the first kind, thus either two conjugate-complex ones or two 
real opposite ones or a real double root, and three roots of the second kind, thus there 


is always one real root. 
Because of the simple character of the moduli of our roots, the procedure of find- 


ing their arguments could be very much abbreviated. However, to show the general 
method, we do not make use of this special property. We begin with the more difficult 
part, namely equation M2. Thus we look back in the chain of transformed equations 
G; until we come to the first equation where our MM, starts to split off. That is Gs, 
where the coefficient of X* is determined to two places. For the transition to G; makes 
two zeros in each of the double products. That is, from Gs an equation M is split off 


1.87X* — 105-4.75X? + 10''-2.08X — 10!7-1.49 = 0 
or, putting X =1.5Y: 
Y* — 0.596Y? + 0.596Y — 1 = (Y — 1)(Y? + 0.404Y + 1). 
This gives 
cos B= — 0.202 or B = 287°.B/32 = 8°.97. 
Thus the argument of the roots equals 
om = 89.97 + 2mr/32, where m = 1, ,32; 


by this the number of values to be tried has sunk from 256 to 32. We can correct 
these values dm by comparing them with the former values $;, that are nearly exact. For 
that we must determine the values of k from the equation ¢,=(2kr+@) /256 = 89.97, 
that is k=~5, thus k=5 and ¢;=8?.9559 - - - 89.956. 

With this value we try to verify the original equation, which on putting x =1.5-y, 


becomes 


16 
or according to (13): 


16 
3 cos 46m + 4 cos 36m + 7 cos 26m + f COS dm + 4 = 0. 


Now the values ¢» are 
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= 89.956, = o1 + 12.59 = 21.4569, --- , = 96°.456 
go = 100° + , bie = 100° + oz, 
giz = 200° + gi, , = 200° + oz, 
= 300% + 1, , = 300% + os. 


These 32 values must be tried. The calculation is carried out to 4-5 decimals. We 
find the solution $3, that is 


= = 1469.4559054397. 


(It is not necessary to compute all 32-4=128 values of cosines, since ¢17, - - - , ds2 
yield values equal or opposite to those yielded by @i, - - - , dis. Also, in the group be- 
longing to ¢; through ¢i¢ not all values are different.) Thus all the roots of the modulus 


1.5 are found. i 
For finding the roots of modulus 1/3, we do not need the somewhat tedious pro- 


cedure above, but can abbreviate it in several ways. A method that is always applica- 
ble consists in dividing the original equation by the product of the three linear factors 


already found, thus by 
(x — 1.5)(x? — 3x cos @ + 2.25). 


In the quotient we put x = y/3 and get y?+1=0, thus y=; so we have as roots of the 
equation 
m= i/3, m= —iV3, = 15, = 
After this general and somewhat tedious way of finding the arguments of the roots 
x4, X5 having the modulus 1.5, we give in the following the simple method appropriate 
to all cases where only two or four complex roots exist. 


From M,, Mz we determine the moduli of their respective roots, as we did above, 
that is X,=X_ and R. Thus the moduli of the roots of the original equation are, as 


above, 
256 


Via td 
and the roots themselves are 
1,2 = V/3-et¥, x3 = X45 = 1.5 -et*, 


Now we use the property of the coefficients of the original equation, namely that 
the sum of all the roots or their combinations at four are: 


+ xo + x3 + t+ x5 = — 0.5 = 20/3 cosy + 1.5 + 3 cos 
— 18/8 = + p/x2+ p/xst+ p/xat p/xs, 


where p = x1%2x%3X4%s5 = 81/8. By inserting the above values we get the two equations 
cos¥ + 3 cos¢@ = — 2, 
2/3 cosy +4 cos = — 8/3, 
cos = — 2/3, or = 146%.4559054397. 


thus 


Then 
cosy = 0 or y = 100%. 
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Conclusions. Graeffe’s method has the following properties. 

i. It yields not only one root at a time but all roots simultaneously, even the com- 
plex ones; this is accomplished by no other method. 

ii. It is the only method that automatically discovers roots lying close together 
which easily escape attention. There is no method other than Graeffe’s which solves, 
without special attention, an equation having, for instance, the roots /5/2 =1.581 -- - 
and 3/2=1.5, not to mention theoretical cases such as 1.67324 --- and 1.67331 --- 
Lagrange’s method is the only other one with the same advantage of separating those 
roots, but the process requires great attention and much computational work and 
fails entirely in the case of complex roots. 

iii. An especially valuable property is that even complex roots of the same modulus 
are automatically obtained. 

iv. These advantages are not due to special artifices. Any other method requires 
a first approximation which must then be corrected. But the finding of this first 
approximation is difficult, particularly in the case of complex roots. Only Bernoulli’s 
method does not require a first approximation, but for that it yields only two real 
roots at a time, and the process of approximation may be very slow. Graeffe’s method 
on the other hand does not require a first approximate value. Besides, it is not neces- 
sary to use criteria of convergence in order to determine if the approximate value is 
sufficiently close the actual root. 

Therefore it seems to us that Graeffe’s method is by far the best for solving algebraic 
equations. Only if one does not need all roots of the equation, but only a single one, 
will it be inferior to other methods. 
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—NOTES— 


ON COMPRESSIBLE FLOW ABOUT BODIES OF REVOLUTION* 
By W. R. SEARS** (Northrop Aircraft, Inc.) 


The linear-perturbation theory of compressible fluid flow originated by Glauert 
and Prandtl has recently been presented in a revised and clarified form by Goldstein 
and Young.' These authors show three alternative procedures by which the compres- 
sible flow in an x, y, z-space can be deduced from a corresponding incompressible flow. 

The linearized differential equation satisfied by the velocity potential @ is 


where 8? denotes 1— U?/a?, U and a being the stream velocity and the velocity of 
sound, respectively, in the undisturbed parallel flow. If the solution of (1) for the 
case 8 =1 (incompressible) is ¢= Ux+f(x, y, z), corresponding solutions for 8B <1 are 
given by the following alternative forms: 


1 
1) @ = Uxt+ By, Bz) 
IT) ¢ = Ux + f(x, By, Bz) 
ITT) = Ux + f(x/B, ¥, 2). 


Each of these variants represents a somewhat different compressible flow, but all 
three are related to the given incompressible flow. The results determined by the the- 
ory are consistent, of course, as far as the linear theory is applicable, and the proced- 
ure used in any given problem is the one that provides the greatest ease of calculation. 
For example, in I the geometry of a slender body remains unaltered as 6 varies; in 
II the body is distorted but the pressures on its surface are unchanged; and so forth. 

Method IT is the one used by Tsien and Lees in a recent paper,’ while both I and 
II are presented by Liepmann and Puckett in a new textbook.’ Sauer‘ writes, in effect, 


IV) @ = Ux + Af(x, By, Bz) 


and selects the value of \ most convenient for any given problem; this includes both I 
and II. Finally, B. Géthert,® rejecting II because of a fancied discrepancy (actually 
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Conclusions. Graeffe’s method has the following properties. 

i. It yields not only one root at a time but all roots simultaneously, even the com- 
plex ones; this is accomplished by no other method. 

ii. It is the only method that automatically discovers roots lying close together 
which easily escape attention. There is no method other than Graeffe’s which solves, 
without special attention, an equation having, for instance, the roots /5/2=1.581 -- - 
and 3/2=1.5, not to mention theoretical cases such as 1.67324 --- and 1.67331 -- - 
Lagrange’s method is the only other one with the same advantage of separating those 
roots, but the process requires great attention and much computational work and 
fails entirely in the case of complex roots. 

iii. An especially valuable property is that even complex roots of the same modulus 
are automatically obtained. 

iv. These advantages are not due to special artifices. Any other method requires 
a first approximation which must then be corrected. But the finding of this first 
approximation is difficult, particularly in the case of complex roots. Only Bernoulli’s 
method does not require a first approximation, but for that it yields only two real 
roots at a time, and the process of approximation may be very slow. Graeffe’s method 
on the other hand does not require a first approximate value. Besides, it is not neces- 
sary to use criteria of convergence in order to determine if the approximate value is 
sufficiently close the actual root. 

Therefore it seems to us that Graeffe’s method is by far the best for solving algebraic 
equations. Only if one does not need all roots of the equation, but only a single one, 
will it be inferior to other methods. 
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and Prandtl has recently been presented in a revised and clarified form by Goldstein 
and Young.! These authors show three alternative procedures by which the compres- 
sible flow in an x, y, z-space can be deduced from a corresponding incompressible flow. 

The linearized differential equation satisfied by the velocity potential ¢ is 


Ox? dy? 02? 
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where 6? denotes 1 — U?/a?, U and a being the stream velocity and the velocity of 
sound, respectively, in the undisturbed parallel flow. If the solution of (1) for the 
case 8=1 (incompressible) is ¢ = Ux+f(x, y, 2), corresponding solutions for 8 <1 are 
given by the following alternative forms: 
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I) Ux + By, Bz) 
II) @ = Ux + f(x, By, Bz) 
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Each of these variants represents a somewhat different compressible flow, but all 
three are related to the given incompressible flow. The results determined by the the- 
ory are consistent, of course, as far as the linear theory is applicable, and the proced- 
ure used in any given problem is the one that provides the greatest ease of calculation. 
For example, in I the geometry of a slender body remains unaltered as B varies; in 
II the body is distorted but the pressures on its surface are unchanged; and so forth. 
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II are presented by Liepmann and Puckett in a new textbook.’ Sauer‘ writes, in effect, 
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5 B. Géthert, Ebene und réumliche Strémung bei hohen Unterschallgeschwindigkeiten, Lilienthal Gesell- 
schaft f. Luftfahrtforschung, Bericht 127, 97-101 (1940). 


191 
4 
sai 
2 ( ) 
— 
4 
— 
a 
= 
2 
= 4 
ts 


192 NOTES [Vol. IV, No. 2 


caused by an error in his application of the method), introduces still another variant 
by writing 
V) @ = BUx + f(x, By, Bz). 


We are particularly concerned here with the application of these procedures to 
the flow about bodies of revolution. There is some confusion on this subject: Gold- 
stein and Young,’ using I, find that the superstream velocities at a streamline body 
in the absence of trailing vortices are 1/8 times those at the same location in incom- 
pressible flow, while both Sauer‘ and Géthert® conclude that these velocities are un- 
affected by compressibility, at least for slender bodies.® 

This confusion is partly due to the fact that, while the procedures are equivalent 
and must yield consistent results in the linear approximation, they may produce dif- 
ferent results when they are applied outside this range. For example, let the maximum 
velocity at the surface of a slender body in incompressible flow be denoted by 
U- [1+ F(n)] where n is the ratio of maximum diameter to length, so that F(n) is a 
given function for any family of bodies. The several variants of the theory then yield 
the following results for the maximum surface velocity (divided by the stream veloc- 


ity) in the compressible flow: 


I) F(n) 
n 
B 
II) F(n/B) 
III) 
IV) 
Vv) Fan). 


Obviously these results are all the same if F(m) is proportional to m or can be 
approximated successfully in that form. But for a typical family of bodies, the ellip- 
soids of revolution, F(m) actually has the form’ 

n? log p — 2n*/1 — n? 1+ — n? 


F(n) = —— where p= — 
— — log p 1—V/1—n? 


(2) 


or, neglecting terms of order n?, 
F(n) = — n? log n. (2a) 
The absence of a linear term in this expression is what leads Géthert to the con- 


clusion that there is no correction for compressibility. Sauer’s similar conclusion 
apparently results somewhat analogously from the fact that he considers only the 


6 Géthert admits a correction “for greater thickness ratio” and proceeds to calculate it by means of 


Method V above. 
7 This is obtained from H. Lamb, Hydrodynamics, Cambridge, 1932, §105. Note that the ratio 


diameter /length, n, is equal to in Lamb's notation. 


| 
| 
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limiting case MO. Actually, in 


(2a) we have retained the leading .20 l 

term while neglecting O(n?), which 

is consistent with the linear-per- @) 

turbation theory. 6 
F(n) according to (2) and (2a) Ln 

is plotted in Figure 1. It is clear o 

that Géthert’s and Sauer’s conclus- 2 

ion cannot be correct in the range of a = v4 

practical interest (1/10<n<1/3 

and 0.6<$6<1, say), since all of -08 Sp 

the various procedures listed above yo 

result in appreciable corrections to os WA 

the velocity ratio. It seems more 


reasonable to conclude merely that 

the linear-perturbation theory can- 0 

not distinguish between the vari- 
0 A 4, 39 

ous results. In this situation the 

formula of Method I might well be Fic. 1. The superstream velocity ratio for ellipsoids 

adopted by reason of its simplicity. of revolution in incompressible flow. 


ON THE NUMERICAL TREATMENT OF FORCED OSCILLATIONS* 
By ALVIN C. SUGAR (Northrop Aircraft) 


1. Introduction. The differential equation, with typical initial conditions, of an 
harmonic oscillator subject to the action of a general disturbing force ma(t) is given by 


wx = a(d), x(0) = 0 = x(0). (1) 


This equation occurs in problems involving from one to infinitely many degrees of 
freedom. Its solution can be expressed as foliows: 


D t 

x =-—», where D -f a(r) sin w(t — r)dr (2) 
0 

is the so-called Duhamel integral. If in (1) we replace only x by D/w we obtain an 

expression for the acceleration of the body. 


# = a(t) — wD. (3) 


In this note a simple expression which is an approximation of D is found. This 
expression provides a convenient process for evaluating x and related quantities. Us- 
ing the resulting simplified form of the acceleration a quick and easy vector method 
of obtaining the maximum acceleration is explained. Rapid methods of finding the 


* Received Oct. 1, 1945. 
** Now at Brown University. 
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maximum displacement are also considered. By maximum acceleration is meant maxi- 
mum magnitude or absolute value of acceleration and similarly for displacement. 
2. Evaluation of the Duhamel integral. The curve y=a(t) is approximated by a 
broken line Y=A(t), see Fig. 1. As indicated we are taking A (to) =0. Since the ap- 
proximation of a curve by a broken line can be improved by increasing the number of 
segments, it is evident that this method can produce a solution which is as accurate 


Y 


t3 


Fic. 1 


as desired. Frequently in engineering problems the forcing function is known only 
approximately and the additional error introduced by a broken line consisting of rela- 
tively few segments is negligible. 

Let A(t) =A,(é), tii StSti, where y;=A,(t) is the equation of a straight line of 
slope Since the lines y;=A,(t) and yi41=Ai4:(¢) intersect at the point A;(t:)), it 
follows that 

= Ai(ti). (4) 
With the notations of Fig. 1 the Duhamel integral can be approximated in the follow- 
ing manner: 


D= J a(r) sin w(t — r)dr 


ty 


A,(r) sin w(t — r)dr + f A ,(r) sin w(t — 7)dr. 


t=] 
Integration by parts enables us to write this in the form 


wD = sal > Mi sin w(t t;_1) > Ki sin w(t 


i=1 i=0 


tj-1 


where we have defined wo=0 in order to obtain the following compact result: 


w i=0 


tm-3 tm 
t_ ths 
tm-2 
t 

| tm-1 

Olt, 
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From (3) and (5) it is evident that 
1H 
A—wD = (wisi — mi) sin w(t — 4). (6) 
t=0 

3. Maximum acceleration and displacement. It is, at times, of importance to know 
the maximum displacement or acceleration. In this paragraph we show a vector 
method of obtaining the maximum acceleration. Since the expression for # for the jth 
time interval is, apart from a constant factor, a sum of sinusoids of the same fre- 
quency, w, it is equal to a sinusoid of frequency w. The amplitude of this sinusoid can 
be obtained by vector methods. It is evident that this amplitude is equal to the maxi- 
mum of the absolute value of the resultant sinusoid over a time interval equal to or 
in excess of one half period. From this it is apparent that the following vector pro- 
cedure can be used in determining max|#| over all of the time intervals. 

If 

— < min ‘= 0, i, (7) 
then max | #| =max OP; (Fig. 2), where the magnitude of the ith vector is 
| (ui—pMis)/w| and its argument with respect to OP; is equal to the phase angle —w;_1. 

Without condition (7) max OP; is an upper bound of | #| and probably a pretty 
good approximation of max | 

Fairly rapid methods of computing the maximum displacement! can be devised 
e.g. when the frequency is large, then, for the 7th interval, max | x| = (1/w?) [max| A,(t)| 
+max| | ]. For any frequency the problem 
of finding max|x| may be reduced to that of P; 
finding the maximum value of the curve ob- 
tained by the superposition of a sinusoid and P, 

a straight line. This can be handled by obvi- sits 
ous methods involving use of elementary = 
differential calculus. ’ 

A still better method of approximating 
the maximum displacement is available if s 
is known or can be quickly evaluated, where s 
is defined by the differential equation, 


s = a(d), $(0) = 0 = s(0). (8) 


This method permits direct use of the above vector procedure. This is easily shown by 
evaluating the Duhamel Integral by repeated integration by parts.? Thus, 


1 t t 
r= —f w(t = f s sin w(t — r)dr 
w Jo 0 


i 
Vx) sin w(t ty), (9) 
kel 
1 A mechanical analyzer was invented by M. A. Biot to obtain this maximum [Bulletin Seismological 


Soc. Amer. 31, 151-171 (1941)]. 
2 As was done by G. W. Housner in obtaining his formulas (2) and (3), Bulletin Seismological Soc. 


Amer, 31, 143-149 (1941). 
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where pv; and ¢, are defined by a broken line approximation of s. If the initial condi- 
tions of (8) defining s are changed to §(0) =so and s(0) =so, then (9) takes the follow- 
ing form: 
So 
r= — sin at socos atts f s sin w(t — r)dr, (10) 
0 
which again can be expressed as a sum of sinusoids of the same frequency. 

On the basis of limited experience, the following suggestions for computation seem 
good. If the curve is sufficiently smooth, then the term containing mw; will make a 
sizeable contribution; consequently the first time interval should be as small as con- 
venient. It seems best to take A(0) =a(0). The vector polygon will obviously be sim- 
plest if t; is selected so that as many values as possible of wt; are multiples of 7. 

If we had assumed A(0) #0, then it would follow that 

wD = A — A(0) cos wt — — x (misa — wi) Sin w(t — 4). 


W 


If we let $,=$(t,) then it is clear that v,41~$,j1. Substituting in (9) we have’ 


2, ($k41 $x) sin w(t ty). (11) 
W k=0 
In calculating the maximum displacement, (11) would be more convenient than (9) 
since § could be obtained by a single integration of a(t). 


A REMARK ON THE RECTIFICATION OF THE 
JOUKOWSKI PROFILE* 


By CHARLES SALTZER (Brown University) 


The Joukowski profile is usually defined as the image under the Joukowski trans- 

formation, 
(1) 

of a circle passing through the point (—c, 0) whose center lies in the first quadrant, 
and whose radius is c(1+¢) where c, and e>0. Although this representation gives the 
complex potential of the incompressible flow about a Joukowski profile very readily, 
the representation of this profile as the inverse of a parabola! has the advantage, as 
will be shown below, of introducing a parameter with direct geometrical meaning 
which permits the immediate rectification of the Joukowski profile in closed form. 

In the 2;-plane consider the parabola 


2 
= (2) 


3 It is interesting to note that the sum in (11) is the so-called left Cauchy-Stieltjes sum corresponding 


to D. 
* Received Aug. 17, 1945. 
1 In this way the profile later called “Joukowski profile” was introduced by Chaplygin. See Chapyl- 


gin’s Collected Papers, Leningrad 1933, vol. 2, pp. 144-178, in particular §6. 
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and the point —(a+07) as a center of inversion. The coefficient of x? may be taken as 
1/2 since it enters only as a scale factor. Only the case a>0, b> —1/2 will be treated 
here.? The transformation 


$1 = (21 + + (3) 


maps the exterior of the parabola in the z;-plane on the exterior of a Joukowski profile 
in the ¢:-plane. For the proper choice of parameters the profiles in the ¢ and {)-planes 
can be mapped on each other by a linear transformation and a reflection. 

Letting ds=|d¢,| and ds;=|dz:| we have for the element of arc length on the 
Joukowski profile by (2) and (3), 

4(1 + x2)!/2 
dx. 

+ + (x? + 26)? 


This expression can be simplified by separation into partial fractions. The roots of 
the denominator can be obtained by equating the latter to zero, transposing one term, 
extracting the square roots of both sides, and solving the two resulting quadratic 
equations. For the non-symmetrical case (a>0) this enables us to write 


= A (= — B, — B,d, 1(< + BB, f, n) |, (5) 


(4) 


ds =|z:+a+ ib|-?| dz, | = 


where (e+ wh 
x 
I ’ ’ ’ = d 6 
(m, p, 9, *1) f (6) 
and 
/2 4a 
a =— [1+ (14+ B=—, = (1+ (7) 
g a 
A = — 86[8* + 48° + a%(6? + d = + + 2)*], f = + — 


The integrand of (6) can be rationalized by the substitution 
x = (1 — u*)/2u (8) 
which gives 


(1 + 2mu — u?)(1 + u?)? 
I(m, p,q, = — 
(m, %1) 2 u?[(u? — pu — 1)? + — p*)u?| 


the same way as the factors of the denominator of (4) were found, and the integration 
can be carried out directly after expanding (9) in partial fractions with linear and 
quadratic denominators. It may be remarked that the J’s taken individually may not 
converge over the entire range of u. 

The case a=0 gives a symmetrical Joukowski profile for which the rectification 
can be carried out in a simpler way. Here equation (4) becomes 


2 This configuration represents the most frequently used Joukowski profiles to within a scale factor 
and a reflection. The case a <0, by reason of symmetry, can be regarded as a reflection of the case a>0, 
and the case bS —1/2 can be treated in a way similar to the treatment of the case b> —1/2. 


? 
} 
= 
where r=4/1+2%1. The factors of the denominator of this integrand can be found in _ 
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2.1/2 


ds = + x3)" + (a1 + 26) | (10) 


and equation (5) is replaced by 


= + a? + (1 — — fa? + (1+ (11) 


Setting 


x= 41 — (12) 
we get, after simplifying and expanding in fractions with quadratic denominators, 
s(x) = g f (13) 
where 
Therefore 


2.—1/2 -1 -1 2,.—1/2 
s(*1) = glo, tanh { (1 + 2x) } —o2 tanh { (1 + x) }]. (15) 
If we denote the slope of the parabola at the point (x1, x{/2) by tan y and note that 
x; =tan y, we can write (15) as 


s(y) = glo. tanh sin y) — o2 (es sin 7) (16) 


where s is measured from the point furthest from the trailing edge. 
In order to introduce the usual parameters ¢ and ¢ for the symmetrical case of the 
Joukowski profile’ consider the circle in the z-plane, 


2(@) = cle + (1+ ee'*] (17) 


which is the image of a symmetrical Joukowski profile in the {-plane. The distance 

between the leading and trailing edges of the Joukowski profile in the {-plane is (re- 

calling Eq. (1)) 

4c(1 + «)? 
1 + 2e 


From (3) the corresponding length in the {1-plane is seen to be 1/6 (i.e. the length 
in the ¢,-plane of the image of the upper half of the imaginary axis in the 2-plane). 
If the profile in the ¢:-plane is identified with the profile in the {-plane then 


b = 3(1 + 26<)(1 + (19) 


If the vertex of the parabola in the a;-plane which corresponds to the given Joukowski 
profile in the ¢-plane is at the origin, and if the x,-axis coincides with the tangent to 
the parabola at this point, then it is readily seen by comparing the positions of the 
profiles in the ¢-plane and the ¢-plane that 


$1 = — if + 2c). (20) 


Since a =0, Eq. (3) can be written 


[2(0)] — ¢[2(x)] = (18) 


1 
ib. (21) 


$1 


3 See, for instance, H. Glauert, The elements of aerofoil and air screw theory, The University Press, 
Cambridge, 1930, pp. 71-75. 
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The successive substitution into Eq. (21) of Eqs. (20), (1), (17), and (19) gives after 
simplification 


tan + i 36. (22) 


1 


€ 
2c(1 + «)? 
This is the parabola 

ce (1 + ©) (23) 


Since a change in the value of c, effects only a change of scale in the {-plane, c may 
be taken without loss of generality as 


c = + (24) 
and this parabola becomes the one considered in Eq. (2). Setting this value of ¢ in 
(19) yields 

b = 4(1 + (25) 
Hence, by Eqs. (7) and (14), 

o=(1 +2) g=/(te. (26) 


Formulae (15) and (16) are valid for o;>0, o2:>0, i.e., for €<1 and thus include those 
profiles whose thicknesses is less than about 4/5 of their lengths. 

It may also be noted that in terms of the variable y, the slope, 0(y), and the 
curvature, d0(y)/ds, for the symmetrical profile may be written as 


4 tan y(tan* y + (1 + 2e)/e*) 
d 
[(1 + 3e)(1 — cos y + 6 cos* y — 3e*/(1 +€)?]. (28) 


CORRECTION AND SUPPLEMENT TO OUR PAPER 

THE CYLINDRICAL ANTENNA: CURRENT AND IMPEDANCE* 
QUARTERLY OF APPLIED MATHEMATICS 3, 302-335 (1946) 
By RONOLD KING anp DAVID MIDDLETON (Harvard University) 


Equation (58) should be written as follows: 


Two lines before this equation | ¥,(0)| /sin Bh should be written instead of | ¥:(0)| , 

These changes involve no alternations in the figures. However, the function 
| ¥(0)| plotted in Fig. 11 to the left of Bk=2/2 is not the parameter of expansion y 
defined by (58) as modified above and as indicated in the caption. The parameter of 
expansion wW as defined in (58) is plotted in Fig. 11a where the part to the right of 
Bh=7/2 is the same as in Fig. 11, the part to the left of 8h=/2 is obtained from the 
curves in Fig. 11 by dividing by sin Bh. 


(58) 


* Received Jan. 25, 1946. 
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For small values of Bh a convenient approximate formula is 
Vxi(0) = 2 — 2 — jph; Bh < 0.5 


so that 


| Wxi(0) | = (2 — 2)? + Bk? = Q — 2 + — 2). 


3.5 


Fic. 11a. The expansion parameter y as defined in the corrected equation (58). 


The following minor errors and misprints have been called to our attention: 

page 312, Eq. (43) change VW to y, 

page 319, Eqs. (59) and 62), change Vx,(z) to W; line following Eq. (61), delete 
the following: y(z) =0 and 

page 320, Eqs. (69) and (70), change b to Q; Eq. (76), insert 1/(m—1)! after the 
first equality sign, ° 

page 323, Eq. (77b), change 4 to y, 

page 324, Eq. (79), insert y after R., 

page 329, Eq. (19), third line, change (Re,;+ 2) to (Re,—w2), 

page 330, Eqs. (23) and (27), page 335 Eqs. (45) and (46), and in the integral 
preceding Eq. (45), change Re, to uz, Ri, to um, throughout, 

page 330, Eq. (24) replace by: u2=(h+2z); u=(h—2), 

page 332, Eq. (43), add superscript bar over first three symbols Ci. 
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